/[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.8 by mlosch, Mon Mar 20 21:36:11 2006 UTC revision 1.16 by mlosch, Tue Apr 24 18:38:15 2007 UTC
# Line 22  C     === Global variables === Line 22  C     === Global variables ===
22  #include "FFIELDS.h"  #include "FFIELDS.h"
23  #include "SEAICE.h"  #include "SEAICE.h"
24  #include "SEAICE_PARAMS.h"  #include "SEAICE_PARAMS.h"
 #include "SEAICE_FFIELDS.h"  
25    
26  C     === Routine arguments ===  C     === Routine arguments ===
27  C     myTime - Simulation time  C     myTime - Simulation time
# Line 39  C     i,j,bi,bj - Loop counters Line 38  C     i,j,bi,bj - Loop counters
38    
39        INTEGER i, j, bi, bj        INTEGER i, j, bi, bj
40        _RL  SINWAT, COSWAT, SINWIN, COSWIN        _RL  SINWAT, COSWAT, SINWIN, COSWIN
41        _RL  fuIce, fvIce, FX, FY        _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)        _RL press       (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
48        _RL etaPlusZeta (1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL sig11       (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
49        _RL zetaMinusEta(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL sig22       (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
50        _RL etaMeanZ    (1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL sig12       (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
51        _RL etaMeanU    (1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL eplus, eminus
       _RL etaMeanV    (1-Olx:sNx+Olx,1-Oly:sNy+Oly)  
       _RL dVdx        (1-Olx:sNx+Olx,1-Oly:sNy+Oly)  
       _RL dVdy        (1-Olx:sNx+Olx,1-Oly:sNy+Oly)  
       _RL dUdx        (1-Olx:sNx+Olx,1-Oly:sNy+Oly)  
       _RL dUdy        (1-Olx:sNx+Olx,1-Oly:sNy+Oly)  
52    
53  c     introduce turning angle (default is zero)  c     introduce turning angle (default is zero)
54        SINWAT=SIN(SEAICE_waterTurnAngle*deg2rad)        SINWAT=SIN(SEAICE_waterTurnAngle*deg2rad)
# Line 82  C Line 79  C
79  C     use an intergral over ice and ocean surface layer to define  C     use an intergral over ice and ocean surface layer to define
80  C     surface stresses on ocean following Hibler and Bryan (1987, JPO)  C     surface stresses on ocean following Hibler and Bryan (1987, JPO)
81  C      C    
82  C     recompute viscosities from updated ice velocities  C     recompute strain rates, viscosities, etc. from updated ice velocities
83         CALL SEAICE_CALC_VISCOSITIES(         IF ( .NOT. SEAICEuseEVP ) THEN
84       I      uIce(1-Olx,1-Oly,1,1,1), vIce(1-Olx,1-Oly,1,1,1),  C     only for EVP we already have the stress components otherwise we need
85       I      zMin, zMax, hEffM, press0,  C     to recompute them here
86       O      eta, zeta, press,          CALL SEAICE_CALC_STRAINRATES(
87  #ifdef SEAICE_ALLOW_EVP       I       uIce(1-Olx,1-Oly,1,1,1), vIce(1-Olx,1-Oly,1,1,1),
88       O      seaice_div, seaice_tension, seaice_shear,       O       e11, e22, e12,
89  #endif /* SEAICE_ALLOW_EVP */       I       myThid )
90       I      myThid )  
91            CALL SEAICE_CALC_VISCOSITIES(
92         I       e11, e22, e12, zMin, zMax, hEffM, press0,
93         O       eta, zeta, press,
94         I       myThid )
95           ENDIF
96  C     re-compute internal stresses with updated ice velocities  C     re-compute internal stresses with updated ice velocities
97         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
98          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
99           DO j=1-Oly+1,sNy+Oly-1           IF ( .NOT. SEAICEuseEVP ) THEN
100            DO i=1-Olx+1,sNx+Olx-1  C     only for EVP we already have computed the stress divergences, for
101             etaPlusZeta (I,J) =  eta(I,J,bi,bj) + zeta(I,J,bi,bj)  C     anything else we have to do it here
102             zetaMinusEta(I,J) = zeta(I,J,bi,bj) -  eta(I,J,bi,bj)            DO j=1-Oly,sNy+Oly
103             etaMeanU (I,J) =             DO i=1-Olx,sNx+Olx
104       &          HALF*(ETA (I,J,bi,bj) + ETA (I-1,J  ,bi,bj))              sig11(I,J) = 0. _d 0
105             etaMeanV (I,J) =              sig22(I,J) = 0. _d 0
106       &          HALF*(ETA (I,J,bi,bj) + ETA (I  ,J-1,bi,bj))              sig12(I,J) = 0. _d 0
107             etaMeanZ (I,J) = QUART *               ENDDO
      &          ( eta(I  ,J,bi,bj) + eta(I  ,J-1,bi,bj)  
      &          + eta(I-1,J,bi,bj) + eta(I-1,J-1,bi,bj) )  
            dUdx(I,J) = ( uIce(I+1,J,1,bi,bj) - uIce(I,J,1,bi,bj) )  
      &          * _recip_dxF(I,J,bi,bj)  
            dUdy(I,J) = ( uIce(I,J+1,1,bi,bj) - uIce(I,J,1,bi,bj) )  
      &          * _recip_dyU(I,J+1,bi,bj)  
            dVdx(I,J) = ( vIce(I+1,J,1,bi,bj) - vIce(I,J,1,bi,bj) )  
      &          * _recip_dxV(I+1,J,bi,bj)  
            dVdy(I,J) = ( vIce(I,J+1,1,bi,bj) - vIce(I,J,1,bi,bj) )  
      &          * _recip_dyF(I,J,bi,bj)  
108            ENDDO            ENDDO
109           ENDDO  
110           DO J = 1,sNy            DO j=1-Oly+1,sNy+Oly-1
111            DO I = 1,sNx             DO i=1-Olx+1,sNx+Olx-1
112  C     First FX = (d/dx)*sigma              eplus = e11(I,J,bi,bj) + e22(I,J,bi,bj)
113  C     + d/dx[ eta+zeta d/dx ] U              eminus= e11(I,J,bi,bj) - e22(I,J,bi,bj)
114             FX = _recip_dxC(I,J,bi,bj) *              sig11(I,J) = zeta(I,J,bi,bj)*eplus + eta(I,J,bi,bj)*eminus
115       &            ( etaPlusZeta(I  ,J) * dUdx(I  ,J)       &           - 0.5 _d 0 * PRESS(I,J,bi,bj)
116       &            - etaPlusZeta(I-1,J) * dUdx(I-1,J) )              sig22(I,J) = zeta(I,J,bi,bj)*eplus - eta(I,J,bi,bj)*eminus
117  C     + (d/dy)[eta*(d/dy + tanphi/a)] U (also on UVRT1/2)       &           - 0.5 _d 0 * PRESS(I,J,bi,bj)
118             FX = FX + _recip_dyG(I,J,bi,bj) * (              sig12(I,J) = 2. _d 0 * e12(I,J,bi,bj) *
119       &          ( etaMeanZ(I,J+1) * dUdy(I,J+1)       &           ( eta(I,J  ,bi,bj) + eta(I-1,J  ,bi,bj)
120       &          - etaMeanZ(I,J  ) * dUdy(I,J  )       &           + eta(I,J-1,bi,bj) + eta(I-1,J-1,bi,bj) )
121       &          )       &           /MAX(1. _d 0,
122       &          - ( etaMeanZ(I,J+1)       &             hEffM(I,J  ,bi,bj) + hEffM(I-1,J  ,bi,bj)
123       &            * ( uIce(I,J+1,1,bi,bj)+uIce(I,J,1,bi,bj) )       &           + hEffM(I,J-1,bi,bj) + hEffM(I-1,J-1,bi,bj))
124       &            - etaMeanZ(I,J  )             ENDDO
125       &            * ( uIce(I,J-1,1,bi,bj)+uIce(I,J,1,bi,bj) ) )            ENDDO
126       &          * 0.5 _d 0 * _tanPhiAtU(I,J,bi,bj)  C     evaluate divergence of stress and apply to forcing
127       &          * recip_rSphere )            DO J=1,sNy
128  C     - 2*eta*(tanphi/a) * ( tanphi/a ) U             DO I=1,sNx
129             FX = FX - TWO * uIce(I,J,1,bi,bj)              FX = ( sig11(I  ,J  ) * _dyF(I  ,J  ,bi,bj)
130       &          * etaMeanU(I,J)*recip_rSphere*recip_rSphere       &           - sig11(I-1,J  ) * _dyF(I-1,J  ,bi,bj)
131       &          * _tanPhiAtU(I,J,bi,bj)  * _tanPhiAtU(I,J,bi,bj)       &           + sig12(I  ,J+1) * _dxV(I  ,J+1,bi,bj)
132  C     + d/dx[ (zeta-eta) dV/dy]       &           - sig12(I  ,J  ) * _dxV(I  ,J  ,bi,bj)
133             FX = FX +       &           ) * recip_rAw(I,J,bi,bj)
134       &          ( zetaMinusEta(I  ,J  ) * dVdy(I  ,J  )       &           -
135       &          - zetaMinusEta(I-1,J  ) * dVdy(I-1,J  )       &           ( sig12(I,J) + sig12(I,J+1) )
136       &          ) * _recip_dxC(I,J,bi,bj)       &           * _tanPhiAtU(I,J,bi,bj) * recip_rSphere
137  C     + d/dy[ eta dV/x ]       &           +
138             FX = FX + (       &           ( sig22(I,J) + sig22(I-1,J) ) * 0.5 _d 0
139       &            etaMeanZ(I,J+1)       &           * _tanPhiAtU(I,J,bi,bj) * recip_rSphere
140       &          * ( vIce(I  ,J+1,1,bi,bj) - vIce(I-1,J+1,1,bi,bj) )  C     one metric term  missing for general curvilinear coordinates
141       &          * _recip_dxV(I,J+1,bi,bj)              FY = ( sig22(I  ,J  ) * _dxF(I  ,J  ,bi,bj)
142       &          - etaMeanZ(I,J  )       &           - sig22(I  ,J-1) * _dxF(I  ,J-1,bi,bj)
143       &          * ( vIce(I  ,J,1,bi,bj) - vIce(I-1,J,1,bi,bj) )       &           + sig12(I+1,J  ) * _dyU(I+1,J  ,bi,bj)
144       &          * _recip_dxV(I,J,bi,bj)       &           - sig12(I  ,J  ) * _dyU(I  ,J  ,bi,bj)
145       &          ) * _recip_dyG(I,J,bi,bj)       &           ) * recip_rAs(I,J,bi,bj)
146  C     - d/dx[ (eta+zeta) * v * (tanphi/a) ]       &           -
147             FX = FX - (       &           ( sig22(I,J) + sig22(I,J-1) ) * 0.5 _d 0
148       &            etaPlusZeta(I  ,J)       &           * _tanPhiAtV(I,J,bi,bj) * recip_rSphere
149       &          * 0.5 * (vIce(I  ,J,1,bi,bj)+vIce(I  ,J+1,1,bi,bj))  C     two metric terms missing for general curvilinear coordinates
      &          * 0.5 * ( _tanPhiAtU(I  ,J,bi,bj)  
      &          + _tanPhiAtU(I+1,J,bi,bj) )  
      &          - etaPlusZeta(I-1,J) *  
      &          * 0.5 * (vIce(I-1,J,1,bi,bj)+vIce(I-1,J+1,1,bi,bj))  
      &          * 0.5 * ( _tanPhiAtU(I-1,J,bi,bj)  
      &          + _tanPhiAtU(I  ,J,bi,bj) )  
      &          )* _recip_dxC(I,J,bi,bj)*recip_rSphere  
 C     - 2*eta*(tanphi/a) * dV/dx  
            FX = FX  
      &          -TWO * etaMeanU(I,J) * _tanPhiAtV(I,J,bi,bj)  
      &          *recip_rSphere  
      &          *(vIce(I  ,J,1,bi,bj) + vIce(I  ,J+1,1,bi,bj)  
      &          - vIce(I-1,J,1,bi,bj) - vIce(I-1,J+1,1,bi,bj))  
      &          * _recip_dxC(I,J,bi,bj)  
 C     - (d/dx) P/2  
            FX = _maskW(I,J,1,bi,bj) * ( FX - _recip_dxC(I,J,bi,bj)  
      &          * ( press(I,J,bi,bj) - press(I-1,J,bi,bj) ) )  
 C  
 C     then FY = (d/dy)*sigma  
 C     + d/dy [(eta+zeta) d/dy] V  
            FY = _recip_dyC(I,J,bi,bj) *  
      &          ( dVdy(I,J  ) * etaPlusZeta(I,J  )  
      &          - dVdy(I,J-1) * etaPlusZeta(I,J-1) )  
 C     + d/dx [eta d/dx] V  
            FY = FY +  _recip_dxC(I,J,bi,bj) *  
      &          ( eta(I  ,J,bi,bj) * dVdx(I  ,J)  
      &          - eta(I-1,J,bi,bj) * dVdx(I-1,J) )  
 C     - d/dy [(zeta-eta) tanphi/a] V  
            FY = FY - _recip_dyC(I,J,bi,bj) * recip_rSphere * (  
      &            zetaMinusEta(I,J  ) * tanPhiAtU(I,J  ,bi,bj)  
      &          * 0.5 * ( vIce(I,J,1,bi,bj) + vIce(I,J+1,1,bi,bj))  
      &          - zetaMinusEta(I,J-1) * tanPhiAtU(I,J-1,bi,bj)  
      &          * 0.5 * ( vIce(I,J,1,bi,bj) + vIce(I,J-1,1,bi,bj)) )  
 C     2*eta tanphi/a ( - tanphi/a - d/dy) V  
            FY = FY - TWO*etaMeanV(I,J) * recip_rSphere  
      &          * _tanPhiAtV(I,J,bi,bj) * (  
      &            _tanPhiAtV(I,J,bi,bj) * recip_rSphere  
      &          + _recip_dyC(I,J,bi,bj) *  
      &          ( 0.5 * ( vIce(I,J,1,bi,bj) + vIce(I,J+1,1,bi,bj))  
      &          - 0.5 * ( vIce(I,J,1,bi,bj) + vIce(I,J-1,1,bi,bj)) ) )  
 C     + d/dy[ (zeta-eta) dU/dx ]  
            FY = FY +  
      &          ( zetaMinusEta(I,J  )*dUdx(I,J  )  
      &          - zetaMinusEta(I,J-1)*dUdx(I,J-1) )  
      &          * _recip_dyC(I,J,bi,bj)  
 C     + d/dx[ eta dU/dy ]  
            FY = FY + _recip_dxG(I,J,bi,bj) *  
      &          ( etaMeanZ(I+1,J) * dUdy(I+1,J)  
      &          - etaMeanZ(I  ,J) * dUdy(I  ,J) )  
 C     + d/dx[ eta * (tanphi/a) * U ]  
            FY = FY + (  
      &            etaMeanZ(I+1,J) * 0.5 *  
      &          ( uIce(I+1,J  ,1,bi,bj) * _tanPhiAtU(I+1,J  ,bi,bj)  
      &          + uIce(I+1,J-1,1,bi,bj) * _tanPhiAtU(I+1,J-1,bi,bj) )  
      &          - etaMeanZ(I  ,J) * 0.5 *  
      &          ( uIce(I  ,J  ,1,bi,bj) * _tanPhiAtU(I  ,J  ,bi,bj)  
      &          + uIce(I  ,J-1,1,bi,bj) * _tanPhiAtU(I  ,J  ,bi,bj) )  
      &          ) *  _recip_dxG(I,J,bi,bj)*recip_rSphere  
 C     + 2*eta*(tanphi/a) dU/dx  
            FY = FY +  
      &          TWO * etaMeanV(I,J)*TWO  * _tanPhiAtV(I,J,bi,bj)  
      &          * ( uIce(I+1,J,1,bi,bj)+uIce(I+1,J-1,1,bi,bj)  
      &            - uIce(I  ,J,1,bi,bj)-uIce(I  ,J-1,1,bi,bj) )  
      &          * _recip_dxG(I,J,bi,bj) * recip_rSphere  
 C     - (d/dy) P/2  
            FY = _maskS(I,J,1,bi,bj) * ( FY - _recip_dyC(I,J,bi,bj)  
      &          * ( press(I,J,bi,bj) - press(I,J-1,bi,bj) ) )  
 C      
 C     recompute wind stress over ice (done already in seaice_dynsolver,  
 C     but not saved)  
            fuIce = 0.5 _d 0 *  
      &          ( DAIRN(I  ,J,bi,bj)*(  
      &          COSWIN*uWind(I  ,J,bi,bj)  
      &          -SIGN(SINWIN, _fCori(I  ,J,bi,bj))*vWind(I  ,J,bi,bj) )  
      &          + DAIRN(I-1,J,bi,bj)*(  
      &          COSWIN*uWind(I-1,J,bi,bj)  
      &          -SIGN(SINWIN, _fCori(I-1,J,bi,bj))*vWind(I-1,J,bi,bj) )  
      &          )  
            fvIce = 0.5 _d 0 *  
      &          ( DAIRN(I,J  ,bi,bj)*(  
      &          SIGN(SINWIN, _fCori(I  ,J,bi,bj))*uWind(I,J  ,bi,bj)  
      &          +COSWIN*vWind(I,J  ,bi,bj) )  
      &          + DAIRN(I,J-1,bi,bj)*(  
      &          SIGN(SINWIN, _fCori(I,J-1,bi,bj))*uWind(I,J-1,bi,bj)  
      &          +COSWIN*vWind(I,J-1,bi,bj) )  
      &          )  
150  C     average wind stress over ice and ocean and apply averaged wind  C     average wind stress over ice and ocean and apply averaged wind
151  C     stress and internal ice stresses to surface layer of ocean  C     stress and internal ice stresses to surface layer of ocean
152             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))
153             areaS = 0.5 * (AREA(I,J,1,bi,bj) + AREA(I,J-1,1,bi,bj))       &           * SEAICEstressFactor
154             fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)+areaW*fuIce + FX              areaS = 0.5 * (AREA(I,J,1,bi,bj) + AREA(I,J-1,1,bi,bj))
155             fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)+areaS*fvIce + FY       &           * SEAICEstressFactor
156            END DO              fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)
157           END DO       &           + areaW*taux(I,J,bi,bj)
158         &           + FX * SEAICEstressFactor
159                fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)
160         &           + areaS*tauy(I,J,bi,bj)
161         &           + FY * SEAICEstressFactor
162    C     save stress divergence for later
163    #ifdef SEAICE_ALLOW_EVP
164                stressDivergenceX(I,J,bi,bj) = FX
165                stressDivergenceY(I,J,bi,bj) = FY
166    #endif /* SEAICE_ALLOW_EVP */
167               ENDDO
168              ENDDO
169             ELSE
170    #ifdef SEAICE_ALLOW_EVP
171              DO J=1,sNy
172               DO I=1,sNx
173    C     average wind stress over ice and ocean and apply averaged wind
174    C     stress and internal ice stresses to surface layer of ocean
175                areaW = 0.5 * (AREA(I,J,1,bi,bj) + AREA(I-1,J,1,bi,bj))
176         &           * SEAICEstressFactor
177                areaS = 0.5 * (AREA(I,J,1,bi,bj) + AREA(I,J-1,1,bi,bj))
178         &           * SEAICEstressFactor
179                fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)
180         &           + areaW*taux(I,J,bi,bj)
181         &           + stressDivergenceX(I,J,bi,bj) * SEAICEstressFactor
182                fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)
183         &           + areaS*tauy(I,J,bi,bj)
184         &           + stressDivergenceY(I,J,bi,bj) * SEAICEstressFactor
185               ENDDO
186              ENDDO
187    #endif /* SEAICE_ALLOW_EVP */
188             ENDIF
189          ENDDO          ENDDO
190         ENDDO         ENDDO
191        ELSE        ELSE
# Line 260  C     ice cover Line 197  C     ice cover
197         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
198          DO j=1,sNy          DO j=1,sNy
199           DO i=1,sNx           DO i=1,sNx
200            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) )*
201       &         COSWAT *       &         COSWAT *
202       &         ( UICE(I,J,1,bi,bj)-GWATX(I,J,bi,bj) )       &         ( UICE(I,J,1,bi,bj)-GWATX(I,J,bi,bj) )
203       &         - SIGN(SINWAT, _fCori(I,J,bi,bj)) * 0.5 _d 0 *       &         - SIGN(SINWAT, _fCori(I,J,bi,bj)) * 0.5 _d 0 *
# Line 271  C     ice cover Line 208  C     ice cover
208       &         0.5 _d 0*(vIce(I-1,J  ,1,bi,bj)-GWATY(I-1,J  ,bi,bj)       &         0.5 _d 0*(vIce(I-1,J  ,1,bi,bj)-GWATY(I-1,J  ,bi,bj)
209       &                  +vIce(I-1,J+1,1,bi,bj)-GWATY(I-1,J+1,bi,bj))       &                  +vIce(I-1,J+1,1,bi,bj)-GWATY(I-1,J+1,bi,bj))
210       &         )       &         )
211            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) )*
212       &         COSWAT *       &         COSWAT *
213       &         ( VICE(I,J,1,bi,bj)-GWATY(I,J,bi,bj) )       &         ( VICE(I,J,1,bi,bj)-GWATY(I,J,bi,bj) )
214       &         + SIGN(SINWAT,  _fCori(I,J,bi,bj)) * 0.5 _d 0 *       &         + SIGN(SINWAT,  _fCori(I,J,bi,bj)) * 0.5 _d 0 *
# Line 283  C     ice cover Line 220  C     ice cover
220       &                  +uIce(I+1,J-1,1,bi,bj)-GWATX(I+1,J-1,bi,bj))       &                  +uIce(I+1,J-1,1,bi,bj)-GWATX(I+1,J-1,bi,bj))
221       &         )       &         )
222            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))
223         &         * SEAICEstressFactor
224            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))
225            fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)+areaW*fuIce       &         * SEAICEstressFactor
226            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
227              fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)+areaS*fvIceLoc
228           ENDDO           ENDDO
229          ENDDO          ENDDO
230         ENDDO         ENDDO
Legend:
Removed from v.1.8  
changed lines
  Added in v.1.16
  ViewVC Help
Powered by ViewVC 1.1.22