/[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.23 by mlosch, Fri May 29 10:18:03 2009 UTC revision 1.25 by mlosch, Wed Jun 24 08:01:43 2009 UTC
# Line 39  C     i,j,bi,bj - Loop counters Line 39  C     i,j,bi,bj - Loop counters
39    
40        INTEGER i, j, bi, bj        INTEGER i, j, bi, bj
41        _RL  SINWAT, COSWAT, SINWIN, COSWIN        _RL  SINWAT, COSWAT, SINWIN, COSWIN
42        _RL  fuIceLoc, fvIceLoc, FX, FY        _RL  fuIceLoc, fvIceLoc
43        _RL  areaW, areaS        _RL  areaW, areaS
44    
       _RL sig11       (1-Olx:sNx+Olx,1-Oly:sNy+Oly)  
       _RL sig22       (1-Olx:sNx+Olx,1-Oly:sNy+Oly)  
       _RL sig12       (1-Olx:sNx+Olx,1-Oly:sNy+Oly)  
       _RL eplus, eminus  
   
45  c     introduce turning angle (default is zero)  c     introduce turning angle (default is zero)
46        SINWAT=SIN(SEAICE_waterTurnAngle*deg2rad)        SINWAT=SIN(SEAICE_waterTurnAngle*deg2rad)
47        COSWAT=COS(SEAICE_waterTurnAngle*deg2rad)        COSWAT=COS(SEAICE_waterTurnAngle*deg2rad)
# Line 58  C Line 53  C
53  C     use an intergral over ice and ocean surface layer to define  C     use an intergral over ice and ocean surface layer to define
54  C     surface stresses on ocean following Hibler and Bryan (1987, JPO)  C     surface stresses on ocean following Hibler and Bryan (1987, JPO)
55  C  C
 C     recompute strain rates, viscosities, etc. from updated ice velocities  
        IF ( .NOT. SEAICEuseEVP ) THEN  
 C     only for EVP we already have the stress components otherwise we need  
 C     to recompute them here  
         CALL SEAICE_CALC_STRAINRATES(  
      I       uIce, vIce,  
      O       e11, e22, e12,  
      I       3, 3, myTime, myIter, myThid )  
   
         CALL SEAICE_CALC_VISCOSITIES(  
      I       e11, e22, e12, zMin, zMax, hEffM, press0,  
      O       eta, zeta, press,  
      I       3, myTime, myIter, myThid )  
        ENDIF  
 C     re-compute internal stresses with updated ice velocities  
56         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
57          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
58           IF ( .NOT. SEAICEuseEVP ) THEN           DO J=1,sNy
59  C     only for EVP we already have computed the stress divergences, for            DO I=1,sNx
 C     anything else we have to do it here  
           DO j=1-Oly,sNy+Oly  
            DO i=1-Olx,sNx+Olx  
             sig11(I,J) = 0. _d 0  
             sig22(I,J) = 0. _d 0  
             sig12(I,J) = 0. _d 0  
            ENDDO  
           ENDDO  
   
           DO j=0,sNy  
            DO i=0,sNx  
             eplus = e11(I,J,bi,bj) + e22(I,J,bi,bj)  
             eminus= e11(I,J,bi,bj) - e22(I,J,bi,bj)  
             sig11(I,J) = zeta(I,J,bi,bj)*eplus + eta(I,J,bi,bj)*eminus  
      &           - 0.5 _d 0 * PRESS(I,J,bi,bj)  
             sig22(I,J) = zeta(I,J,bi,bj)*eplus - eta(I,J,bi,bj)*eminus  
      &           - 0.5 _d 0 * PRESS(I,J,bi,bj)  
            ENDDO  
           ENDDO  
   
           DO j=1,sNy+1  
            DO i=1,sNx+1  
             sig12(I,J) = 2. _d 0 * e12(I,J,bi,bj) *  
      &           ( eta(I,J  ,bi,bj) + eta(I-1,J  ,bi,bj)  
      &           + eta(I,J-1,bi,bj) + eta(I-1,J-1,bi,bj) )  
      &           /MAX(1. _d 0,  
      &             hEffM(I,J  ,bi,bj) + hEffM(I-1,J  ,bi,bj)  
      &           + hEffM(I,J-1,bi,bj) + hEffM(I-1,J-1,bi,bj))  
            ENDDO  
           ENDDO  
 C     evaluate divergence of stress and apply to forcing  
           DO J=1,sNy  
            DO I=1,sNx  
             FX = ( sig11(I  ,J  ) * _dyF(I  ,J  ,bi,bj)  
      &           - sig11(I-1,J  ) * _dyF(I-1,J  ,bi,bj)  
      &           + sig12(I  ,J+1) * _dxV(I  ,J+1,bi,bj)  
      &           - sig12(I  ,J  ) * _dxV(I  ,J  ,bi,bj)  
      &           ) * recip_rAw(I,J,bi,bj)  
             FY = ( sig22(I  ,J  ) * _dxF(I  ,J  ,bi,bj)  
      &           - sig22(I  ,J-1) * _dxF(I  ,J-1,bi,bj)  
      &           + sig12(I+1,J  ) * _dyU(I+1,J  ,bi,bj)  
      &           - sig12(I  ,J  ) * _dyU(I  ,J  ,bi,bj)  
      &           ) * recip_rAs(I,J,bi,bj)  
 C     average wind stress over ice and ocean and apply averaged wind  
 C     stress and internal ice stresses to surface layer of ocean  
             areaW = 0.5 * (AREA(I,J,1,bi,bj) + AREA(I-1,J,1,bi,bj))  
      &           * SEAICEstressFactor  
             areaS = 0.5 * (AREA(I,J,1,bi,bj) + AREA(I,J-1,1,bi,bj))  
      &           * SEAICEstressFactor  
             fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)  
      &           + areaW*taux(I,J,bi,bj)  
      &           + FX * SEAICEstressFactor  
             fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)  
      &           + areaS*tauy(I,J,bi,bj)  
      &           + FY * SEAICEstressFactor  
 C     save stress divergence for later  
 #ifdef SEAICE_ALLOW_EVP  
             stressDivergenceX(I,J,bi,bj) = FX  
             stressDivergenceY(I,J,bi,bj) = FY  
 #endif /* SEAICE_ALLOW_EVP */  
            ENDDO  
           ENDDO  
          ELSE  
 #ifdef SEAICE_ALLOW_EVP  
           DO J=1,sNy  
            DO I=1,sNx  
60  C     average wind stress over ice and ocean and apply averaged wind  C     average wind stress over ice and ocean and apply averaged wind
61  C     stress and internal ice stresses to surface layer of ocean  C     stress and internal ice stresses to surface layer of ocean
62              areaW = 0.5 * (AREA(I,J,1,bi,bj) + AREA(I-1,J,1,bi,bj))             areaW = 0.5 * (AREA(I,J,bi,bj) + AREA(I-1,J,bi,bj))
63       &           * SEAICEstressFactor       &          * SEAICEstressFactor
64              areaS = 0.5 * (AREA(I,J,1,bi,bj) + AREA(I,J-1,1,bi,bj))             areaS = 0.5 * (AREA(I,J,bi,bj) + AREA(I,J-1,bi,bj))
65       &           * SEAICEstressFactor       &          * SEAICEstressFactor
66              fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)             fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)
67       &           + areaW*taux(I,J,bi,bj)       &          + areaW*taux(I,J,bi,bj)
68       &           + stressDivergenceX(I,J,bi,bj) * SEAICEstressFactor       &          + stressDivergenceX(I,J,bi,bj) * SEAICEstressFactor
69              fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)             fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)
70       &           + areaS*tauy(I,J,bi,bj)       &          + areaS*tauy(I,J,bi,bj)
71       &           + stressDivergenceY(I,J,bi,bj) * SEAICEstressFactor       &          + stressDivergenceY(I,J,bi,bj) * SEAICEstressFactor
            ENDDO  
72            ENDDO            ENDDO
73  #endif /* SEAICE_ALLOW_EVP */           ENDDO
          ENDIF  
74          ENDDO          ENDDO
75         ENDDO         ENDDO
76    
# Line 174  C     ice cover Line 86  C     ice cover
86           DO i=1,sNx           DO i=1,sNx
87            fuIceLoc=HALF*( DWATN(I,J,bi,bj)+DWATN(I-1,J,bi,bj) )*            fuIceLoc=HALF*( DWATN(I,J,bi,bj)+DWATN(I-1,J,bi,bj) )*
88       &         COSWAT *       &         COSWAT *
89       &         ( UICE(I,J,1,bi,bj)-uVel(I,J,1,bi,bj) )       &         ( UICE(I,J,bi,bj)-uVel(I,J,1,bi,bj) )
90       &         - SIGN(SINWAT, _fCori(I,J,bi,bj)) * 0.5 _d 0 *       &         - SIGN(SINWAT, _fCori(I,J,bi,bj)) * 0.5 _d 0 *
91       &         ( DWATN(I  ,J,bi,bj) *       &         ( DWATN(I  ,J,bi,bj) *
92       &         0.5 _d 0*(vIce(I  ,J  ,1,bi,bj)-vVel(I  ,J  ,1,bi,bj)       &         0.5 _d 0*(vIce(I  ,J  ,bi,bj)-vVel(I  ,J  ,1,bi,bj)
93       &                  +vIce(I  ,J+1,1,bi,bj)-vVel(I  ,J+1,1,bi,bj))       &                  +vIce(I  ,J+1,bi,bj)-vVel(I  ,J+1,1,bi,bj))
94       &         + DWATN(I-1,J,bi,bj) *       &         + DWATN(I-1,J,bi,bj) *
95       &         0.5 _d 0*(vIce(I-1,J  ,1,bi,bj)-vVel(I-1,J  ,1,bi,bj)       &         0.5 _d 0*(vIce(I-1,J  ,bi,bj)-vVel(I-1,J  ,1,bi,bj)
96       &                  +vIce(I-1,J+1,1,bi,bj)-vVel(I-1,J+1,1,bi,bj))       &                  +vIce(I-1,J+1,bi,bj)-vVel(I-1,J+1,1,bi,bj))
97       &         )       &         )
98            fvIceLoc=HALF*( DWATN(I,J,bi,bj)+DWATN(I,J-1,bi,bj) )*            fvIceLoc=HALF*( DWATN(I,J,bi,bj)+DWATN(I,J-1,bi,bj) )*
99       &         COSWAT *       &         COSWAT *
100       &         ( VICE(I,J,1,bi,bj)-vVel(I,J,1,bi,bj) )       &         ( VICE(I,J,bi,bj)-vVel(I,J,1,bi,bj) )
101       &         + SIGN(SINWAT,  _fCori(I,J,bi,bj)) * 0.5 _d 0 *       &         + SIGN(SINWAT,  _fCori(I,J,bi,bj)) * 0.5 _d 0 *
102       &         ( DWATN(I,J  ,bi,bj) *       &         ( DWATN(I,J  ,bi,bj) *
103       &         0.5 _d 0*(uIce(I  ,J  ,1,bi,bj)-uVel(I  ,J  ,1,bi,bj)       &         0.5 _d 0*(uIce(I  ,J  ,bi,bj)-uVel(I  ,J  ,1,bi,bj)
104       &                  +uIce(I+1,J  ,1,bi,bj)-uVel(I+1,J  ,1,bi,bj))       &                  +uIce(I+1,J  ,bi,bj)-uVel(I+1,J  ,1,bi,bj))
105       &         + DWATN(I,J-1,bi,bj) *       &         + DWATN(I,J-1,bi,bj) *
106       &         0.5 _d 0*(uIce(I  ,J-1,1,bi,bj)-uVel(I  ,J-1,1,bi,bj)       &         0.5 _d 0*(uIce(I  ,J-1,bi,bj)-uVel(I  ,J-1,1,bi,bj)
107       &                  +uIce(I+1,J-1,1,bi,bj)-uVel(I+1,J-1,1,bi,bj))       &                  +uIce(I+1,J-1,bi,bj)-uVel(I+1,J-1,1,bi,bj))
108       &         )       &         )
109            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,bi,bj) + AREA(I-1,J,bi,bj))
110       &         * SEAICEstressFactor       &         * SEAICEstressFactor
111            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,bi,bj) + AREA(I,J-1,bi,bj))
112       &         * SEAICEstressFactor       &         * SEAICEstressFactor
113            fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)+areaW*fuIceLoc            fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)+areaW*fuIceLoc
114            fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)+areaS*fvIceLoc            fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)+areaS*fvIceLoc
Legend:
Removed from v.1.23  
changed lines
  Added in v.1.25
  ViewVC Help
Powered by ViewVC 1.1.22