/[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.18 by mlosch, Tue May 15 14:32:56 2007 UTC revision 1.19 by jmc, Tue Nov 13 19:26:26 2007 UTC
# Line 4  C $Name$ Line 4  C $Name$
4  #include "SEAICE_OPTIONS.h"  #include "SEAICE_OPTIONS.h"
5    
6  CStartOfInterface  CStartOfInterface
7        SUBROUTINE SEAICE_OCEAN_STRESS(        SUBROUTINE SEAICE_OCEAN_STRESS(
8       I     myTime, myIter, myThid )       I     myTime, myIter, myThid )
9  C     /==========================================================\  C     /==========================================================\
10  C     | SUBROUTINE SEAICE_OCEAN_STRESS                           |  C     | SUBROUTINE SEAICE_OCEAN_STRESS                           |
# Line 32  C     myThid - Thread no. that called th Line 32  C     myThid - Thread no. that called th
32        INTEGER myThid        INTEGER myThid
33  CEndOfInterface  CEndOfInterface
34    
35  #ifdef SEAICE_CGRID  #ifdef SEAICE_CGRID
36  C     === Local variables ===  C     === Local variables ===
37  C     i,j,bi,bj - Loop counters  C     i,j,bi,bj - Loop counters
38    
# Line 58  c     introduce turning angle (default i Line 58  c     introduce turning angle (default i
58    
59        IF ( useHB87StressCoupling ) THEN        IF ( useHB87StressCoupling ) THEN
60  C  C
61  C     use an intergral over ice and ocean surface layer to define  C     use an intergral over ice and ocean surface layer to define
62  C     surface stresses on ocean following Hibler and Bryan (1987, JPO)  C     surface stresses on ocean following Hibler and Bryan (1987, JPO)
63  C      C
64  C     recompute strain rates, viscosities, etc. from updated ice velocities  C     recompute strain rates, viscosities, etc. from updated ice velocities
65         IF ( .NOT. SEAICEuseEVP ) THEN         IF ( .NOT. SEAICEuseEVP ) THEN
66  C     only for EVP we already have the stress components otherwise we need  C     only for EVP we already have the stress components otherwise we need
67  C     to recompute them here  C     to recompute them here
68          CALL SEAICE_CALC_STRAINRATES(          CALL SEAICE_CALC_STRAINRATES(
69       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),
70       O       e11, e22, e12,       O       e11, e22, e12,
71       I       myThid )       I       3, myTime, myIter, myThid )
72    
73          CALL SEAICE_CALC_VISCOSITIES(          CALL SEAICE_CALC_VISCOSITIES(
74       I       e11, e22, e12, zMin, zMax, hEffM, press0,       I       e11, e22, e12, zMin, zMax, hEffM, press0,
75       O       eta, zeta, press,       O       eta, zeta, press,
76       I       myThid )       I       3, myTime, myIter, myThid )
77         ENDIF         ENDIF
78  C     re-compute internal stresses with updated ice velocities  C     re-compute internal stresses with updated ice velocities
79         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
80          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
81           IF ( .NOT. SEAICEuseEVP ) THEN           IF ( .NOT. SEAICEuseEVP ) THEN
82  C     only for EVP we already have computed the stress divergences, for  C     only for EVP we already have computed the stress divergences, for
83  C     anything else we have to do it here  C     anything else we have to do it here
84            DO j=1-Oly,sNy+Oly            DO j=1-Oly,sNy+Oly
85             DO i=1-Olx,sNx+Olx             DO i=1-Olx,sNx+Olx
# Line 105  C     anything else we have to do it her Line 105  C     anything else we have to do it her
105              sig12(I,J) = 2. _d 0 * e12(I,J,bi,bj) *              sig12(I,J) = 2. _d 0 * e12(I,J,bi,bj) *
106       &           ( eta(I,J  ,bi,bj) + eta(I-1,J  ,bi,bj)       &           ( eta(I,J  ,bi,bj) + eta(I-1,J  ,bi,bj)
107       &           + eta(I,J-1,bi,bj) + eta(I-1,J-1,bi,bj) )       &           + eta(I,J-1,bi,bj) + eta(I-1,J-1,bi,bj) )
108       &           /MAX(1. _d 0,       &           /MAX(1. _d 0,
109       &             hEffM(I,J  ,bi,bj) + hEffM(I-1,J  ,bi,bj)       &             hEffM(I,J  ,bi,bj) + hEffM(I-1,J  ,bi,bj)
110       &           + hEffM(I,J-1,bi,bj) + hEffM(I-1,J-1,bi,bj))       &           + hEffM(I,J-1,bi,bj) + hEffM(I-1,J-1,bi,bj))
111             ENDDO             ENDDO
# Line 115  C     evaluate divergence of stress and Line 115  C     evaluate divergence of stress and
115             DO I=1,sNx             DO I=1,sNx
116              FX = ( sig11(I  ,J  ) * _dyF(I  ,J  ,bi,bj)              FX = ( sig11(I  ,J  ) * _dyF(I  ,J  ,bi,bj)
117       &           - sig11(I-1,J  ) * _dyF(I-1,J  ,bi,bj)       &           - sig11(I-1,J  ) * _dyF(I-1,J  ,bi,bj)
118       &           + sig12(I  ,J+1) * _dxV(I  ,J+1,bi,bj)       &           + sig12(I  ,J+1) * _dxV(I  ,J+1,bi,bj)
119       &           - sig12(I  ,J  ) * _dxV(I  ,J  ,bi,bj)       &           - sig12(I  ,J  ) * _dxV(I  ,J  ,bi,bj)
120       &           ) * recip_rAw(I,J,bi,bj)       &           ) * recip_rAw(I,J,bi,bj)
121       &           -       &           -
122       &           ( sig12(I,J) + sig12(I,J+1) )       &           ( sig12(I,J) + sig12(I,J+1) )
123       &           * _tanPhiAtU(I,J,bi,bj) * recip_rSphere       &           * _tanPhiAtU(I,J,bi,bj) * recip_rSphere
124       &           +       &           +
125       &           ( sig22(I,J) + sig22(I-1,J) ) * 0.5 _d 0       &           ( sig22(I,J) + sig22(I-1,J) ) * 0.5 _d 0
126       &           * _tanPhiAtU(I,J,bi,bj) * recip_rSphere       &           * _tanPhiAtU(I,J,bi,bj) * recip_rSphere
127  C     one metric term  missing for general curvilinear coordinates  C     one metric term  missing for general curvilinear coordinates
128              FY = ( sig22(I  ,J  ) * _dxF(I  ,J  ,bi,bj)              FY = ( sig22(I  ,J  ) * _dxF(I  ,J  ,bi,bj)
129       &           - sig22(I  ,J-1) * _dxF(I  ,J-1,bi,bj)       &           - sig22(I  ,J-1) * _dxF(I  ,J-1,bi,bj)
130       &           + sig12(I+1,J  ) * _dyU(I+1,J  ,bi,bj)       &           + sig12(I+1,J  ) * _dyU(I+1,J  ,bi,bj)
131       &           - sig12(I  ,J  ) * _dyU(I  ,J  ,bi,bj)       &           - sig12(I  ,J  ) * _dyU(I  ,J  ,bi,bj)
132       &           ) * recip_rAs(I,J,bi,bj)       &           ) * recip_rAs(I,J,bi,bj)
133       &           -       &           -
134       &           ( sig22(I,J) + sig22(I,J-1) ) * 0.5 _d 0       &           ( sig22(I,J) + sig22(I,J-1) ) * 0.5 _d 0
135       &           * _tanPhiAtV(I,J,bi,bj) * recip_rSphere       &           * _tanPhiAtV(I,J,bi,bj) * recip_rSphere
136  C     two metric terms missing for general curvilinear coordinates  C     two metric terms missing for general curvilinear coordinates
137  C     average wind stress over ice and ocean and apply averaged wind  C     average wind stress over ice and ocean and apply averaged wind
138  C     stress and internal ice stresses to surface layer of ocean  C     stress and internal ice stresses to surface layer of ocean
139              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))
140       &           * SEAICEstressFactor       &           * SEAICEstressFactor
141              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))
142       &           * SEAICEstressFactor       &           * SEAICEstressFactor
143              fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)              fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)
144       &           + areaW*taux(I,J,bi,bj)       &           + areaW*taux(I,J,bi,bj)
145       &           + FX * SEAICEstressFactor       &           + FX * SEAICEstressFactor
146              fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)              fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)
147       &           + areaS*tauy(I,J,bi,bj)       &           + areaS*tauy(I,J,bi,bj)
# Line 157  C     save stress divergence for later Line 157  C     save stress divergence for later
157  #ifdef SEAICE_ALLOW_EVP  #ifdef SEAICE_ALLOW_EVP
158            DO J=1,sNy            DO J=1,sNy
159             DO I=1,sNx             DO I=1,sNx
160  C     average wind stress over ice and ocean and apply averaged wind  C     average wind stress over ice and ocean and apply averaged wind
161  C     stress and internal ice stresses to surface layer of ocean  C     stress and internal ice stresses to surface layer of ocean
162              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))
163       &           * SEAICEstressFactor       &           * SEAICEstressFactor
164              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))
165       &           * SEAICEstressFactor       &           * SEAICEstressFactor
166              fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)              fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)
167       &           + areaW*taux(I,J,bi,bj)       &           + areaW*taux(I,J,bi,bj)
168       &           + stressDivergenceX(I,J,bi,bj) * SEAICEstressFactor       &           + stressDivergenceX(I,J,bi,bj) * SEAICEstressFactor
169              fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)              fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)
170       &           + areaS*tauy(I,J,bi,bj)       &           + areaS*tauy(I,J,bi,bj)
# Line 179  C     stress and internal ice stresses t Line 179  C     stress and internal ice stresses t
179        ELSE        ELSE
180  C     else: useHB87StressCoupling=F  C     else: useHB87StressCoupling=F
181    
182  C--   Compute ice-affected wind stress (interpolate to U/V-points)  C--   Compute ice-affected wind stress (interpolate to U/V-points)
183  C     by averaging wind stress and ice-ocean stress according to  C     by averaging wind stress and ice-ocean stress according to
184  C     ice cover  C     ice cover
185        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
186         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
187          DO j=1,sNy          DO j=1,sNy
188           DO i=1,sNx           DO i=1,sNx
189            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) )*
190       &         COSWAT *       &         COSWAT *
191       &         ( UICE(I,J,1,bi,bj)-GWATX(I,J,bi,bj) )       &         ( UICE(I,J,1,bi,bj)-GWATX(I,J,bi,bj) )
192       &         - SIGN(SINWAT, _fCori(I,J,bi,bj)) * 0.5 _d 0 *       &         - SIGN(SINWAT, _fCori(I,J,bi,bj)) * 0.5 _d 0 *
193       &         ( DWATN(I  ,J,bi,bj) *       &         ( DWATN(I  ,J,bi,bj) *
194       &         0.5 _d 0*(vIce(I  ,J  ,1,bi,bj)-GWATY(I  ,J  ,bi,bj)       &         0.5 _d 0*(vIce(I  ,J  ,1,bi,bj)-GWATY(I  ,J  ,bi,bj)
195       &                  +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))
196       &         + DWATN(I-1,J,bi,bj) *       &         + DWATN(I-1,J,bi,bj) *
197       &         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)
198       &                  +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))
199       &         )       &         )
200            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) )*
201       &         COSWAT *       &         COSWAT *
# Line 205  C     ice cover Line 205  C     ice cover
205       &         0.5 _d 0*(uIce(I  ,J  ,1,bi,bj)-GWATX(I  ,J  ,bi,bj)       &         0.5 _d 0*(uIce(I  ,J  ,1,bi,bj)-GWATX(I  ,J  ,bi,bj)
206       &                  +uIce(I+1,J  ,1,bi,bj)-GWATX(I+1,J  ,bi,bj))       &                  +uIce(I+1,J  ,1,bi,bj)-GWATX(I+1,J  ,bi,bj))
207       &         + DWATN(I,J-1,bi,bj) *       &         + DWATN(I,J-1,bi,bj) *
208       &         0.5 _d 0*(uIce(I  ,J-1,1,bi,bj)-GWATX(I  ,J-1,bi,bj)       &         0.5 _d 0*(uIce(I  ,J-1,1,bi,bj)-GWATX(I  ,J-1,bi,bj)
209       &                  +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))
210       &         )       &         )
211            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))
212       &         * SEAICEstressFactor       &         * SEAICEstressFactor
Legend:
Removed from v.1.18  
changed lines
  Added in v.1.19
  ViewVC Help
Powered by ViewVC 1.1.22