/[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.5 by mlosch, Wed Mar 15 19:49:04 2006 UTC revision 1.8 by mlosch, Mon Mar 20 21:36:11 2006 UTC
# Line 87  C     recompute viscosities from updated Line 87  C     recompute viscosities from updated
87       I      uIce(1-Olx,1-Oly,1,1,1), vIce(1-Olx,1-Oly,1,1,1),       I      uIce(1-Olx,1-Oly,1,1,1), vIce(1-Olx,1-Oly,1,1,1),
88       I      zMin, zMax, hEffM, press0,       I      zMin, zMax, hEffM, press0,
89       O      eta, zeta, press,       O      eta, zeta, press,
90    #ifdef SEAICE_ALLOW_EVP
91         O      seaice_div, seaice_tension, seaice_shear,
92    #endif /* SEAICE_ALLOW_EVP */
93       I      myThid )       I      myThid )
94  C     re-compute internal stresses with updated ice velocities  C     re-compute internal stresses with updated ice velocities
95         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
# Line 114  C     re-compute internal stresses with Line 117  C     re-compute internal stresses with
117           ENDDO           ENDDO
118           DO J = 1,sNy           DO J = 1,sNy
119            DO I = 1,sNx            DO I = 1,sNx
120  C     First FX = (d/dx)*sigam  C     First FX = (d/dx)*sigma
121  C     + d/dx[ eta+zeta d/dx ] U  C     + d/dx[ eta+zeta d/dx ] U
122             FX = _recip_dxC(I,J,bi,bj) *             FX = _recip_dxC(I,J,bi,bj) *
123       &            ( etaPlusZeta(I  ,J) * dUdx(I  ,J)       &            ( etaPlusZeta(I  ,J) * dUdx(I  ,J)
# Line 170  C     - (d/dx) P/2 Line 173  C     - (d/dx) P/2
173             FX = _maskW(I,J,1,bi,bj) * ( FX - _recip_dxC(I,J,bi,bj)             FX = _maskW(I,J,1,bi,bj) * ( FX - _recip_dxC(I,J,bi,bj)
174       &          * ( press(I,J,bi,bj) - press(I-1,J,bi,bj) ) )       &          * ( press(I,J,bi,bj) - press(I-1,J,bi,bj) ) )
175  C  C
176  C     then FY = (d/dy)*sigam  C     then FY = (d/dy)*sigma
177  C     + d/dy [(eta+zeta) d/dy] V  C     + d/dy [(eta+zeta) d/dy] V
178             FY = _recip_dyC(I,J,bi,bj) *             FY = _recip_dyC(I,J,bi,bj) *
179       &          ( dVdy(I,J  ) * etaPlusZeta(I,J  )       &          ( dVdy(I,J  ) * etaPlusZeta(I,J  )
# Line 224  C     recompute wind stress over ice (do Line 227  C     recompute wind stress over ice (do
227  C     but not saved)  C     but not saved)
228             fuIce = 0.5 _d 0 *             fuIce = 0.5 _d 0 *
229       &          ( DAIRN(I  ,J,bi,bj)*(       &          ( DAIRN(I  ,J,bi,bj)*(
230       &          COSWIN*uWind(I  ,J,bi,bj)-SINWIN*vWind(I  ,J,bi,bj) )       &          COSWIN*uWind(I  ,J,bi,bj)
231         &          -SIGN(SINWIN, _fCori(I  ,J,bi,bj))*vWind(I  ,J,bi,bj) )
232       &          + DAIRN(I-1,J,bi,bj)*(       &          + DAIRN(I-1,J,bi,bj)*(
233       &          COSWIN*uWind(I-1,J,bi,bj)-SINWIN*vWind(I-1,J,bi,bj) )       &          COSWIN*uWind(I-1,J,bi,bj)
234         &          -SIGN(SINWIN, _fCori(I-1,J,bi,bj))*vWind(I-1,J,bi,bj) )
235       &          )       &          )
236             fvIce = 0.5 _d 0 *             fvIce = 0.5 _d 0 *
237       &          ( DAIRN(I,J  ,bi,bj)*(       &          ( DAIRN(I,J  ,bi,bj)*(
238       &          SINWIN*uWind(I,J  ,bi,bj)+COSWIN*vWind(I,J  ,bi,bj) )       &          SIGN(SINWIN, _fCori(I  ,J,bi,bj))*uWind(I,J  ,bi,bj)
239         &          +COSWIN*vWind(I,J  ,bi,bj) )
240       &          + DAIRN(I,J-1,bi,bj)*(       &          + DAIRN(I,J-1,bi,bj)*(
241       &          SINWIN*uWind(I,J-1,bi,bj)+COSWIN*vWind(I,J-1,bi,bj) )       &          SIGN(SINWIN, _fCori(I,J-1,bi,bj))*uWind(I,J-1,bi,bj)
242         &          +COSWIN*vWind(I,J-1,bi,bj) )
243       &          )       &          )
244  C     average wind stress over ice and ocean and apply averaged wind  C     average wind stress over ice and ocean and apply averaged wind
245  C     stress and internal ice stresses to surface layer of ocean  C     stress and internal ice stresses to surface layer of ocean
246             areaW = 0.5 * (AREA(I,J,1,bi,bj) + AREA(I-1,J,1,bi,bj))             areaW = 0.5 * (AREA(I,J,1,bi,bj) + AREA(I-1,J,1,bi,bj))
247             areaS = 0.5 * (AREA(I,J,1,bi,bj) + AREA(I,J-1,1,bi,bj))             areaS = 0.5 * (AREA(I,J,1,bi,bj) + AREA(I,J-1,1,bi,bj))
248             fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)+areaW*(fuIce + FX)             fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)+areaW*fuIce + FX
249             fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)+areaS*(fvIce + FY)             fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)+areaS*fvIce + FY
250            END DO            END DO
251           END DO           END DO
252          ENDDO          ENDDO
# Line 253  C     ice cover Line 260  C     ice cover
260         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
261          DO j=1,sNy          DO j=1,sNy
262           DO i=1,sNx           DO i=1,sNx
263            fuIce=HALF*( DWATN(I,J,bi,bj)+DWATN(I,J+1,bi,bj) )*(            fuIce=HALF*( DWATN(I,J,bi,bj)+DWATN(I,J+1,bi,bj) )*
264       &         COSWAT *       &         COSWAT *
265       &         ( UICE(I,J,1,bi,bj)-GWATX(I,J,bi,bj) )       &         ( UICE(I,J,1,bi,bj)-GWATX(I,J,bi,bj) )
266       &         - SINWAT* 0.5 _d 0 * (       &         - SIGN(SINWAT, _fCori(I,J,bi,bj)) * 0.5 _d 0 *
267       &          0.5 _d 0*(vIce(I  ,J  ,1,bi,bj)-GWATY(I  ,J  ,bi,bj)       &         ( DWATN(I  ,J,bi,bj) *
268       &                   +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)
269       &         +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))
270       &                   +vIce(I-1,J+1,1,bi,bj)-GWATY(I-1,J+1,bi,bj)) )       &         + DWATN(I-1,J,bi,bj) *
271         &         0.5 _d 0*(vIce(I-1,J  ,1,bi,bj)-GWATY(I-1,J  ,bi,bj)
272         &                  +vIce(I-1,J+1,1,bi,bj)-GWATY(I-1,J+1,bi,bj))
273       &         )       &         )
274            fvIce=HALF*( DWATN(I,J,bi,bj)+DWATN(I+1,J,bi,bj) )*(            fvIce=HALF*( DWATN(I,J,bi,bj)+DWATN(I+1,J,bi,bj) )*
275       &         SINWAT *       &         COSWAT *
276       &         ( UICE(I,J,1,bi,bj)-GWATX(I,J,bi,bj) )       &         ( VICE(I,J,1,bi,bj)-GWATY(I,J,bi,bj) )
277       &         + COSWAT * 0.5 _d 0 * (       &         + SIGN(SINWAT,  _fCori(I,J,bi,bj)) * 0.5 _d 0 *
278       &          0.5 _d 0*(uIce(I  ,J  ,1,bi,bj)-GWATX(I  ,J  ,bi,bj)       &         ( DWATN(I,J  ,bi,bj) *
279       &                   +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)
280       &         +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))
281       &                   +uIce(I+1,J-1,1,bi,bj)-GWATX(I+1,J-1,bi,bj)) )       &         + DWATN(I,J-1,bi,bj) *
282         &         0.5 _d 0*(uIce(I  ,J-1,1,bi,bj)-GWATX(I  ,J-1,bi,bj)
283         &                  +uIce(I+1,J-1,1,bi,bj)-GWATX(I+1,J-1,bi,bj))
284       &         )       &         )
285            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))
286            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))
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.22