/[MITgcm]/MITgcm/pkg/seaice/seaice_ocean_stress.F
ViewVC logotype

Diff of /MITgcm/pkg/seaice/seaice_ocean_stress.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.4 by mlosch, Thu Mar 9 20:22:40 2006 UTC revision 1.13 by mlosch, Fri Apr 20 18:29:58 2007 UTC
# Line 18  C     === Global variables === Line 18  C     === Global variables ===
18  #include "SIZE.h"  #include "SIZE.h"
19  #include "EEPARAMS.h"  #include "EEPARAMS.h"
20  #include "PARAMS.h"  #include "PARAMS.h"
21    #include "GRID.h"
22  #include "FFIELDS.h"  #include "FFIELDS.h"
23  #include "SEAICE.h"  #include "SEAICE.h"
24  #include "SEAICE_PARAMS.h"  #include "SEAICE_PARAMS.h"
# Line 29  C     myThid - Thread no. that called th Line 30  C     myThid - Thread no. that called th
30        _RL     myTime        _RL     myTime
31        INTEGER myIter        INTEGER myIter
32        INTEGER myThid        INTEGER myThid
 CML      _RL COR_ICE    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,  nSx,nSy)  
33  CEndOfInterface  CEndOfInterface
34    
35  #ifdef SEAICE_CGRID  #ifdef SEAICE_CGRID
# Line 37  C     === Local variables === Line 37  C     === Local variables ===
37  C     i,j,bi,bj - Loop counters  C     i,j,bi,bj - Loop counters
38    
39        INTEGER i, j, bi, bj        INTEGER i, j, bi, bj
40        _RL  SINWAT, COSWAT        _RL  SINWAT, COSWAT, SINWIN, COSWIN
41        _RL  fuIce, fvIce        _RL  fuIceLoc, fvIceLoc, FX, FY
42        _RL  areaW, areaS        _RL  areaW, areaS
43    
44          _RL e11         (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
45          _RL e22         (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
46          _RL e12         (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
47          _RL press       (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
48          _RL etaPlusZeta (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
49          _RL zetaMinusEta(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
50          _RL etaMeanZ    (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
51          _RL etaMeanU    (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
52          _RL etaMeanV    (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
53          _RL dVdx        (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
54          _RL dVdy        (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
55          _RL dUdx        (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
56          _RL dUdy        (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
57    
58  c     introduce turning angle (default is zero)  c     introduce turning angle (default is zero)
59        SINWAT=SIN(SEAICE_waterTurnAngle*deg2rad)        SINWAT=SIN(SEAICE_waterTurnAngle*deg2rad)
60        COSWAT=COS(SEAICE_waterTurnAngle*deg2rad)        COSWAT=COS(SEAICE_waterTurnAngle*deg2rad)
61          SINWIN=SIN(SEAICE_airTurnAngle*deg2rad)
62          COSWIN=COS(SEAICE_airTurnAngle*deg2rad)
63    
64  C--   Update overlap regions  C--   Update overlap regions
65        CALL EXCH_UV_XY_RL(WINDX, WINDY, .TRUE., myThid)        CALL EXCH_UV_XY_RL(WINDX, WINDY, .TRUE., myThid)
# Line 63  C     to U and V points of C-grid for fo Line 79  C     to U and V points of C-grid for fo
79        ENDDO        ENDDO
80  #endif /* ifndef SEAICE_EXTERNAL_FLUXES */  #endif /* ifndef SEAICE_EXTERNAL_FLUXES */
81    
82  C--   Compute ice-affected wind stress (interpolate to U/V-points)        IF ( useHB87StressCoupling ) THEN
83    C
84    C     use an intergral over ice and ocean surface layer to define
85    C     surface stresses on ocean following Hibler and Bryan (1987, JPO)
86    C    
87    C     recompute viscosities from updated ice velocities
88          CALL SEAICE_CALC_STRAINRATES(
89         I     uIce(1-Olx,1-Oly,1,1,1), vIce(1-Olx,1-Oly,1,1,1),
90         O     e11, e22, e12,
91         I     myThid )
92    
93          CALL SEAICE_CALC_VISCOSITIES(
94         I     e11, e22, e12, zMin, zMax, hEffM, press0,
95         O     eta, zeta, press,
96         I     myThid )
97    C     re-compute internal stresses with updated ice velocities
98           DO bj=myByLo(myThid),myByHi(myThid)
99            DO bi=myBxLo(myThid),myBxHi(myThid)
100             DO j=1-Oly+1,sNy+Oly-1
101              DO i=1-Olx+1,sNx+Olx-1
102               etaPlusZeta (I,J) =  eta(I,J,bi,bj) + zeta(I,J,bi,bj)
103               zetaMinusEta(I,J) = zeta(I,J,bi,bj) -  eta(I,J,bi,bj)
104               etaMeanU (I,J) =
105         &          HALF*(ETA (I,J,bi,bj) + ETA (I-1,J  ,bi,bj))
106               etaMeanV (I,J) =
107         &          HALF*(ETA (I,J,bi,bj) + ETA (I  ,J-1,bi,bj))
108               etaMeanZ (I,J) = QUART *  
109         &          ( eta(I  ,J,bi,bj) + eta(I  ,J-1,bi,bj)
110         &          + eta(I-1,J,bi,bj) + eta(I-1,J-1,bi,bj) )
111               dUdx(I,J) = ( uIce(I+1,J,1,bi,bj) - uIce(I,J,1,bi,bj) )
112         &          * _recip_dxF(I,J,bi,bj)
113               dUdy(I,J) = ( uIce(I,J+1,1,bi,bj) - uIce(I,J,1,bi,bj) )
114         &          * _recip_dyU(I,J+1,bi,bj)
115               dVdx(I,J) = ( vIce(I+1,J,1,bi,bj) - vIce(I,J,1,bi,bj) )
116         &          * _recip_dxV(I+1,J,bi,bj)
117               dVdy(I,J) = ( vIce(I,J+1,1,bi,bj) - vIce(I,J,1,bi,bj) )
118         &          * _recip_dyF(I,J,bi,bj)
119              ENDDO
120             ENDDO
121             DO J = 1,sNy
122              DO I = 1,sNx
123    C     First FX = (d/dx)*sigma
124    C     + d/dx[ eta+zeta d/dx ] U
125               FX = _recip_dxC(I,J,bi,bj) *
126         &            ( etaPlusZeta(I  ,J) * dUdx(I  ,J)
127         &            - etaPlusZeta(I-1,J) * dUdx(I-1,J) )
128    C     + (d/dy)[eta*(d/dy + tanphi/a)] U (also on UVRT1/2)
129               FX = FX + _recip_dyG(I,J,bi,bj) * (
130         &          ( etaMeanZ(I,J+1) * dUdy(I,J+1)
131         &          - etaMeanZ(I,J  ) * dUdy(I,J  )
132         &          )
133         &          - ( etaMeanZ(I,J+1)
134         &            * ( uIce(I,J+1,1,bi,bj)+uIce(I,J,1,bi,bj) )
135         &            - etaMeanZ(I,J  )
136         &            * ( uIce(I,J-1,1,bi,bj)+uIce(I,J,1,bi,bj) ) )
137         &          * 0.5 _d 0 * _tanPhiAtU(I,J,bi,bj)
138         &          * recip_rSphere )
139    C     - 2*eta*(tanphi/a) * ( tanphi/a ) U
140               FX = FX - TWO * uIce(I,J,1,bi,bj)
141         &          * etaMeanU(I,J)*recip_rSphere*recip_rSphere
142         &          * _tanPhiAtU(I,J,bi,bj)  * _tanPhiAtU(I,J,bi,bj)
143    C     + d/dx[ (zeta-eta) dV/dy]
144               FX = FX +
145         &          ( zetaMinusEta(I  ,J  ) * dVdy(I  ,J  )
146         &          - zetaMinusEta(I-1,J  ) * dVdy(I-1,J  )
147         &          ) * _recip_dxC(I,J,bi,bj)
148    C     + d/dy[ eta dV/x ]
149               FX = FX + (
150         &            etaMeanZ(I,J+1)
151         &          * ( vIce(I  ,J+1,1,bi,bj) - vIce(I-1,J+1,1,bi,bj) )
152         &          * _recip_dxV(I,J+1,bi,bj)
153         &          - etaMeanZ(I,J  )
154         &          * ( vIce(I  ,J,1,bi,bj) - vIce(I-1,J,1,bi,bj) )
155         &          * _recip_dxV(I,J,bi,bj)
156         &          ) * _recip_dyG(I,J,bi,bj)
157    C     - d/dx[ (eta+zeta) * v * (tanphi/a) ]
158               FX = FX - (
159         &            etaPlusZeta(I  ,J)
160         &          * 0.5 * (vIce(I  ,J,1,bi,bj)+vIce(I  ,J+1,1,bi,bj))
161         &          * 0.5 * ( _tanPhiAtU(I  ,J,bi,bj)
162         &          + _tanPhiAtU(I+1,J,bi,bj) )
163         &          - etaPlusZeta(I-1,J) *
164         &          * 0.5 * (vIce(I-1,J,1,bi,bj)+vIce(I-1,J+1,1,bi,bj))
165         &          * 0.5 * ( _tanPhiAtU(I-1,J,bi,bj)
166         &          + _tanPhiAtU(I  ,J,bi,bj) )
167         &          )* _recip_dxC(I,J,bi,bj)*recip_rSphere
168    C     - 2*eta*(tanphi/a) * dV/dx
169               FX = FX
170         &          -TWO * etaMeanU(I,J) * _tanPhiAtV(I,J,bi,bj)
171         &          *recip_rSphere
172         &          *(vIce(I  ,J,1,bi,bj) + vIce(I  ,J+1,1,bi,bj)
173         &          - vIce(I-1,J,1,bi,bj) - vIce(I-1,J+1,1,bi,bj))
174         &          * _recip_dxC(I,J,bi,bj)
175    C     - (d/dx) P/2
176               FX = _maskW(I,J,1,bi,bj) * ( FX - _recip_dxC(I,J,bi,bj)
177         &          * ( press(I,J,bi,bj) - press(I-1,J,bi,bj) ) )
178    C
179    C     then FY = (d/dy)*sigma
180    C     + d/dy [(eta+zeta) d/dy] V
181               FY = _recip_dyC(I,J,bi,bj) *
182         &          ( dVdy(I,J  ) * etaPlusZeta(I,J  )
183         &          - dVdy(I,J-1) * etaPlusZeta(I,J-1) )
184    C     + d/dx [eta d/dx] V
185               FY = FY +  _recip_dxC(I,J,bi,bj) *
186         &          ( eta(I  ,J,bi,bj) * dVdx(I  ,J)
187         &          - eta(I-1,J,bi,bj) * dVdx(I-1,J) )
188    C     - d/dy [(zeta-eta) tanphi/a] V
189               FY = FY - _recip_dyC(I,J,bi,bj) * recip_rSphere * (
190         &            zetaMinusEta(I,J  ) * tanPhiAtU(I,J  ,bi,bj)
191         &          * 0.5 * ( vIce(I,J,1,bi,bj) + vIce(I,J+1,1,bi,bj))
192         &          - zetaMinusEta(I,J-1) * tanPhiAtU(I,J-1,bi,bj)
193         &          * 0.5 * ( vIce(I,J,1,bi,bj) + vIce(I,J-1,1,bi,bj)) )
194    C     2*eta tanphi/a ( - tanphi/a - d/dy) V
195               FY = FY - TWO*etaMeanV(I,J) * recip_rSphere
196         &          * _tanPhiAtV(I,J,bi,bj) * (
197         &            _tanPhiAtV(I,J,bi,bj) * recip_rSphere
198         &          + _recip_dyC(I,J,bi,bj) *
199         &          ( 0.5 * ( vIce(I,J,1,bi,bj) + vIce(I,J+1,1,bi,bj))
200         &          - 0.5 * ( vIce(I,J,1,bi,bj) + vIce(I,J-1,1,bi,bj)) ) )
201    C     + d/dy[ (zeta-eta) dU/dx ]
202               FY = FY +
203         &          ( zetaMinusEta(I,J  )*dUdx(I,J  )
204         &          - zetaMinusEta(I,J-1)*dUdx(I,J-1) )
205         &          * _recip_dyC(I,J,bi,bj)
206    C     + d/dx[ eta dU/dy ]
207               FY = FY + _recip_dxG(I,J,bi,bj) *
208         &          ( etaMeanZ(I+1,J) * dUdy(I+1,J)
209         &          - etaMeanZ(I  ,J) * dUdy(I  ,J) )
210    C     + d/dx[ eta * (tanphi/a) * U ]
211               FY = FY + (
212         &            etaMeanZ(I+1,J) * 0.5 *
213         &          ( uIce(I+1,J  ,1,bi,bj) * _tanPhiAtU(I+1,J  ,bi,bj)
214         &          + uIce(I+1,J-1,1,bi,bj) * _tanPhiAtU(I+1,J-1,bi,bj) )
215         &          - etaMeanZ(I  ,J) * 0.5 *
216         &          ( uIce(I  ,J  ,1,bi,bj) * _tanPhiAtU(I  ,J  ,bi,bj)
217         &          + uIce(I  ,J-1,1,bi,bj) * _tanPhiAtU(I  ,J  ,bi,bj) )
218         &          ) *  _recip_dxG(I,J,bi,bj)*recip_rSphere
219    C     + 2*eta*(tanphi/a) dU/dx
220               FY = FY +
221         &          TWO * etaMeanV(I,J)*TWO  * _tanPhiAtV(I,J,bi,bj)
222         &          * ( uIce(I+1,J,1,bi,bj)+uIce(I+1,J-1,1,bi,bj)
223         &            - uIce(I  ,J,1,bi,bj)-uIce(I  ,J-1,1,bi,bj) )
224         &          * _recip_dxG(I,J,bi,bj) * recip_rSphere
225    C     - (d/dy) P/2
226               FY = _maskS(I,J,1,bi,bj) * ( FY - _recip_dyC(I,J,bi,bj)
227         &          * ( press(I,J,bi,bj) - press(I,J-1,bi,bj) ) )
228    C     average wind stress over ice and ocean and apply averaged wind
229    C     stress and internal ice stresses to surface layer of ocean
230               areaW = 0.5 * (AREA(I,J,1,bi,bj) + AREA(I-1,J,1,bi,bj))
231         &         * SEAICEstressFactor
232               areaS = 0.5 * (AREA(I,J,1,bi,bj) + AREA(I,J-1,1,bi,bj))
233         &         * SEAICEstressFactor
234               fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)
235         &          + areaW*taux(I,J,bi,bj)
236         &          + FX * SEAICEstressFactor
237               fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)
238         &          + areaS*tauy(I,J,bi,bj)
239         &          + FY * SEAICEstressFactor
240              END DO
241             END DO
242            ENDDO
243           ENDDO
244          ELSE
245    
246    C--   Compute ice-affected wind stress (interpolate to U/V-points)
247    C     by averaging wind stress and ice-ocean stress according to
248    C     ice cover
249        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
250         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
251          DO j=1,sNy          DO j=1,sNy
252           DO i=1,sNx           DO i=1,sNx
253            fuIce=HALF*( DWATN(I,J,bi,bj)+DWATN(I,J+1,bi,bj) )*(            fuIceLoc=HALF*( DWATN(I,J,bi,bj)+DWATN(I,J+1,bi,bj) )*
254       &         COSWAT *       &         COSWAT *
255       &         ( UICE(I,J,1,bi,bj)-GWATX(I,J,bi,bj) )       &         ( UICE(I,J,1,bi,bj)-GWATX(I,J,bi,bj) )
256       &         - SINWAT* 0.5 _d 0 * (       &         - SIGN(SINWAT, _fCori(I,J,bi,bj)) * 0.5 _d 0 *
257       &          0.5 _d 0*(vIce(I  ,J  ,1,bi,bj)-GWATY(I  ,J  ,bi,bj)       &         ( DWATN(I  ,J,bi,bj) *
258       &                   +vIce(I-1,J  ,1,bi,bj)-GWATY(I-1,J  ,bi,bj))       &         0.5 _d 0*(vIce(I  ,J  ,1,bi,bj)-GWATY(I  ,J  ,bi,bj)
259       &         +0.5 _d 0*(vIce(I  ,J+1,1,bi,bj)-GWATY(I  ,J+1,bi,bj)       &                  +vIce(I  ,J+1,1,bi,bj)-GWATY(I  ,J+1,bi,bj))
260       &                   +vIce(I-1,J+1,1,bi,bj)-GWATY(I-1,J+1,bi,bj)) )       &         + DWATN(I-1,J,bi,bj) *
261         &         0.5 _d 0*(vIce(I-1,J  ,1,bi,bj)-GWATY(I-1,J  ,bi,bj)
262         &                  +vIce(I-1,J+1,1,bi,bj)-GWATY(I-1,J+1,bi,bj))
263       &         )       &         )
264            fvIce=HALF*( DWATN(I,J,bi,bj)+DWATN(I+1,J,bi,bj) )*(            fvIceLoc=HALF*( DWATN(I,J,bi,bj)+DWATN(I+1,J,bi,bj) )*
265       &         SINWAT *       &         COSWAT *
266       &         ( UICE(I,J,1,bi,bj)-GWATX(I,J,bi,bj) )       &         ( VICE(I,J,1,bi,bj)-GWATY(I,J,bi,bj) )
267       &         + COSWAT * 0.5 _d 0 * (       &         + SIGN(SINWAT,  _fCori(I,J,bi,bj)) * 0.5 _d 0 *
268       &          0.5 _d 0*(uIce(I  ,J  ,1,bi,bj)-GWATX(I  ,J  ,bi,bj)       &         ( DWATN(I,J  ,bi,bj) *
269       &                   +uIce(I+1,J  ,1,bi,bj)-GWATX(I+1,J  ,bi,bj))       &         0.5 _d 0*(uIce(I  ,J  ,1,bi,bj)-GWATX(I  ,J  ,bi,bj)
270       &         +0.5 _d 0*(uIce(I  ,J-1,1,bi,bj)-GWATX(I  ,J-1,bi,bj)       &                  +uIce(I+1,J  ,1,bi,bj)-GWATX(I+1,J  ,bi,bj))
271       &                   +uIce(I+1,J-1,1,bi,bj)-GWATX(I+1,J-1,bi,bj)) )       &         + DWATN(I,J-1,bi,bj) *
272         &         0.5 _d 0*(uIce(I  ,J-1,1,bi,bj)-GWATX(I  ,J-1,bi,bj)
273         &                  +uIce(I+1,J-1,1,bi,bj)-GWATX(I+1,J-1,bi,bj))
274       &         )       &         )
275            areaW = 0.5 _d 0 * (AREA(I,J,1,bi,bj) + AREA(I-1,J,1,bi,bj))            areaW = 0.5 _d 0 * (AREA(I,J,1,bi,bj) + AREA(I-1,J,1,bi,bj))
276         &         * SEAICEstressFactor
277            areaS = 0.5 _d 0 * (AREA(I,J,1,bi,bj) + AREA(I,J-1,1,bi,bj))            areaS = 0.5 _d 0 * (AREA(I,J,1,bi,bj) + AREA(I,J-1,1,bi,bj))
278            fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)+areaW*fuIce       &         * SEAICEstressFactor
279            fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)+areaS*fvIce            fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)+areaW*fuIceLoc
280              fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)+areaS*fvIceLoc
281           ENDDO           ENDDO
282          ENDDO          ENDDO
283         ENDDO         ENDDO
284        ENDDO        ENDDO
285          ENDIF
286        CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)        CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)
287    
288  #endif /* not SEAICE_CGRID */  #endif /* not SEAICE_CGRID */
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.13
  ViewVC Help
Powered by ViewVC 1.1.22