/[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.1 by mlosch, Mon Mar 6 13:17:37 2006 UTC revision 1.21 by dimitri, Thu Jan 17 23:18:39 2008 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 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 "DYNVARS.h"
22    #include "GRID.h"
23  #include "FFIELDS.h"  #include "FFIELDS.h"
24  #include "SEAICE.h"  #include "SEAICE.h"
25  #include "SEAICE_PARAMS.h"  #include "SEAICE_PARAMS.h"
# Line 29  C     myThid - Thread no. that called th Line 31  C     myThid - Thread no. that called th
31        _RL     myTime        _RL     myTime
32        INTEGER myIter        INTEGER myIter
33        INTEGER myThid        INTEGER myThid
 CML      _RL COR_ICE    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,  nSx,nSy)  
34  CEndOfInterface  CEndOfInterface
35    
36  #ifdef SEAICE_CGRID  #ifdef SEAICE_CGRID
37  C     === Local variables ===  C     === Local variables ===
38  C     i,j,bi,bj - Loop counters  C     i,j,bi,bj - Loop counters
39    
40        INTEGER i, j, bi, bj        INTEGER i, j, bi, bj
41        _RL  SINWIN, COSWIN, SINWAT, COSWAT        _RL  SINWAT, COSWAT, SINWIN, COSWIN
42  #ifdef SEAICE_TEST_ICE_STRESS_1        _RL  fuIceLoc, fvIceLoc, FX, FY
43        _RL  fuIce, fvIce        _RL  areaW, areaS
44  #endif  
45          _RL e11         (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
46          _RL e22         (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
47          _RL e12         (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
48          _RL press       (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
49          _RL sig11       (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
50          _RL sig22       (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
51          _RL sig12       (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
52          _RL eplus, eminus
53    
54  c     introduce turning angle (default is zero)  c     introduce turning angle (default is zero)
       SINWIN=SIN(SEAICE_airTurnAngle*deg2rad)  
       COSWIN=COS(SEAICE_airTurnAngle*deg2rad)  
55        SINWAT=SIN(SEAICE_waterTurnAngle*deg2rad)        SINWAT=SIN(SEAICE_waterTurnAngle*deg2rad)
56        COSWAT=COS(SEAICE_waterTurnAngle*deg2rad)        COSWAT=COS(SEAICE_waterTurnAngle*deg2rad)
57          SINWIN=SIN(SEAICE_airTurnAngle*deg2rad)
58          COSWIN=COS(SEAICE_airTurnAngle*deg2rad)
59    
60  CML#ifdef SEAICE_ORIGINAL_BAD_ICE_STRESS        IF ( useHB87StressCoupling ) THEN
61  CMLC--   Following formulation is problematic and is no longer used.  C
62  CML#ifdef SEAICE_ALLOW_DYNAMICS  C     use an intergral over ice and ocean surface layer to define
63  CML      IF ( SEAICEuseDYNAMICS ) THEN  C     surface stresses on ocean following Hibler and Bryan (1987, JPO)
64  CMLC--   Compute ice-affected wind stress  C
65  CML       DO bj=myByLo(myThid),myByHi(myThid)  C     recompute strain rates, viscosities, etc. from updated ice velocities
66  CML        DO bi=myBxLo(myThid),myBxHi(myThid)         IF ( .NOT. SEAICEuseEVP ) THEN
67  CML         DO j=1,sNy  C     only for EVP we already have the stress components otherwise we need
68  CML          DO i=1,sNx  C     to recompute them here
69  CML           WINDX(I,J,bi,bj)=DWATN(I,J,bi,bj)          CALL SEAICE_CALC_STRAINRATES(
70  CML     &          *(COSWAT*(GWATX(I,J,bi,bj)-UICE(I,J,1,bi,bj))       I       uIce, vIce,
71  CML     &          -SINWAT*(GWATY(I,J,bi,bj)-VICEC(I,J,bi,bj)))       O       e11, e22, e12,
72  CML           WINDY(I,J,bi,bj)=DWATN(I,J,bi,bj)       I       3, 3, myTime, myIter, myThid )
73  CML     &          *(SINWAT*(GWATX(I,J,bi,bj)-UICEC(I,J,bi,bj))  
74  CML     &          +COSWAT*(GWATY(I,J,bi,bj)-VICE(I,J,1,bi,bj)))          CALL SEAICE_CALC_VISCOSITIES(
75  CML           WINDX(I,J,bi,bj)=WINDX(I,J,bi,bj)-( COR_ICE(I,J,bi,bj)       I       e11, e22, e12, zMin, zMax, hEffM, press0,
76  CML     &          *GWATY(I,J,bi,bj)-COR_ICE(I,J,bi,bj)*VICEC(I,J,bi,bj))       O       eta, zeta, press,
77  CML           WINDY(I,J,bi,bj)=WINDY(I,J,bi,bj)-(-COR_ICE(I,J,bi,bj)       I       3, myTime, myIter, myThid )
78  CML     &          *GWATX(I,J,bi,bj)+COR_ICE(I,J,bi,bj)*UICEC(I,J,bi,bj))         ENDIF
79  CML           WINDX(I,J,bi,bj)=WINDX(I,J,bi,bj)-(UICE(I,J,1,bi,bj)  C     re-compute internal stresses with updated ice velocities
80  CML     &          -UICE(I,J,3,bi,bj))*AMASS(I,J,bi,bj)/SEAICE_DT*TWO         DO bj=myByLo(myThid),myByHi(myThid)
81  CML           WINDY(I,J,bi,bj)=WINDY(I,J,bi,bj)-(VICE(I,J,1,bi,bj)          DO bi=myBxLo(myThid),myBxHi(myThid)
82  CML     &          -VICE(I,J,3,bi,bj))*AMASS(I,J,bi,bj)/SEAICE_DT*TWO           IF ( .NOT. SEAICEuseEVP ) THEN
83  CML          ENDDO  C     only for EVP we already have computed the stress divergences, for
84  CML         ENDDO  C     anything else we have to do it here
85  CML        ENDDO            DO j=1-Oly,sNy+Oly
86  CML       ENDDO             DO i=1-Olx,sNx+Olx
87  CML       DO bj=myByLo(myThid),myByHi(myThid)              sig11(I,J) = 0. _d 0
88  CML        DO bi=myBxLo(myThid),myBxHi(myThid)              sig22(I,J) = 0. _d 0
89  CML         DO j=1,sNy              sig12(I,J) = 0. _d 0
90  CML          DO i=1,sNx             ENDDO
91  CML           WINDX(I,J,bi,bj)=-WINDX(I,J,bi,bj)            ENDDO
92  CML           WINDY(I,J,bi,bj)=-WINDY(I,J,bi,bj)  
93  CML          ENDDO            DO j=0,sNy
94  CML         ENDDO             DO i=0,sNx
95  CML        ENDDO              eplus = e11(I,J,bi,bj) + e22(I,J,bi,bj)
96  CML       ENDDO              eminus= e11(I,J,bi,bj) - e22(I,J,bi,bj)
97  CML      ENDIF              sig11(I,J) = zeta(I,J,bi,bj)*eplus + eta(I,J,bi,bj)*eminus
98  CML#endif /*  SEAICE_ALLOW_DYNAMICS */       &           - 0.5 _d 0 * PRESS(I,J,bi,bj)
99  CML#endif /* SEAICE_ORIGINAL_BAD_ICE_STRESS */              sig22(I,J) = zeta(I,J,bi,bj)*eplus - eta(I,J,bi,bj)*eminus
100         &           - 0.5 _d 0 * PRESS(I,J,bi,bj)
101  C--   Update overlap regions             ENDDO
102        CALL EXCH_UV_XY_RL(WINDX, WINDY, .TRUE., myThid)            ENDDO
103    
104  #ifndef SEAICE_EXTERNAL_FLUXES            DO j=1,sNy+1
105  C--   Interpolate wind stress (N/m^2) from South-West B-grid             DO i=1,sNx+1
106  C     to South-West C-grid for forcing ocean model.              sig12(I,J) = 2. _d 0 * e12(I,J,bi,bj) *
107        DO bj=myByLo(myThid),myByHi(myThid)       &           ( eta(I,J  ,bi,bj) + eta(I-1,J  ,bi,bj)
108         DO bi=myBxLo(myThid),myBxHi(myThid)       &           + eta(I,J-1,bi,bj) + eta(I-1,J-1,bi,bj) )
109          DO j=1,sNy       &           /MAX(1. _d 0,
110           DO i=1,sNx       &             hEffM(I,J  ,bi,bj) + hEffM(I-1,J  ,bi,bj)
111              fu(I,J,bi,bj)=WINDX(I,J,bi,bj)       &           + hEffM(I,J-1,bi,bj) + hEffM(I-1,J-1,bi,bj))
112              fv(I,J,bi,bj)=WINDY(I,J,bi,bj)             ENDDO
113           ENDDO            ENDDO
114    C     evaluate divergence of stress and apply to forcing
115              DO J=1,sNy
116               DO I=1,sNx
117                FX = ( sig11(I  ,J  ) * _dyF(I  ,J  ,bi,bj)
118         &           - sig11(I-1,J  ) * _dyF(I-1,J  ,bi,bj)
119         &           + sig12(I  ,J+1) * _dxV(I  ,J+1,bi,bj)
120         &           - sig12(I  ,J  ) * _dxV(I  ,J  ,bi,bj)
121         &           ) * recip_rAw(I,J,bi,bj)
122         &           -
123         &           ( sig12(I,J) + sig12(I,J+1) )
124         &           * _tanPhiAtU(I,J,bi,bj) * recip_rSphere
125         &           +
126         &           ( sig22(I,J) + sig22(I-1,J) ) * 0.5 _d 0
127         &           * _tanPhiAtU(I,J,bi,bj) * recip_rSphere
128    C     one metric term  missing for general curvilinear coordinates
129                FY = ( sig22(I  ,J  ) * _dxF(I  ,J  ,bi,bj)
130         &           - sig22(I  ,J-1) * _dxF(I  ,J-1,bi,bj)
131         &           + sig12(I+1,J  ) * _dyU(I+1,J  ,bi,bj)
132         &           - sig12(I  ,J  ) * _dyU(I  ,J  ,bi,bj)
133         &           ) * recip_rAs(I,J,bi,bj)
134         &           -
135         &           ( sig22(I,J) + sig22(I,J-1) ) * 0.5 _d 0
136         &           * _tanPhiAtV(I,J,bi,bj) * recip_rSphere
137    C     two metric terms missing for general curvilinear coordinates
138    C     average wind stress over ice and ocean and apply averaged wind
139    C     stress and internal ice stresses to surface layer of ocean
140                areaW = 0.5 * (AREA(I,J,1,bi,bj) + AREA(I-1,J,1,bi,bj))
141         &           * SEAICEstressFactor
142                areaS = 0.5 * (AREA(I,J,1,bi,bj) + AREA(I,J-1,1,bi,bj))
143         &           * SEAICEstressFactor
144                fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)
145         &           + areaW*taux(I,J,bi,bj)
146         &           + FX * SEAICEstressFactor
147                fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)
148         &           + areaS*tauy(I,J,bi,bj)
149         &           + FY * SEAICEstressFactor
150    C     save stress divergence for later
151    #ifdef SEAICE_ALLOW_EVP
152                stressDivergenceX(I,J,bi,bj) = FX
153                stressDivergenceY(I,J,bi,bj) = FY
154    #endif /* SEAICE_ALLOW_EVP */
155               ENDDO
156              ENDDO
157             ELSE
158    #ifdef SEAICE_ALLOW_EVP
159              DO J=1,sNy
160               DO I=1,sNx
161    C     average wind stress over ice and ocean and apply averaged wind
162    C     stress and internal ice stresses to surface layer of ocean
163                areaW = 0.5 * (AREA(I,J,1,bi,bj) + AREA(I-1,J,1,bi,bj))
164         &           * SEAICEstressFactor
165                areaS = 0.5 * (AREA(I,J,1,bi,bj) + AREA(I,J-1,1,bi,bj))
166         &           * SEAICEstressFactor
167                fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)
168         &           + areaW*taux(I,J,bi,bj)
169         &           + stressDivergenceX(I,J,bi,bj) * SEAICEstressFactor
170                fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)
171         &           + areaS*tauy(I,J,bi,bj)
172         &           + stressDivergenceY(I,J,bi,bj) * SEAICEstressFactor
173               ENDDO
174              ENDDO
175    #endif /* SEAICE_ALLOW_EVP */
176             ENDIF
177          ENDDO          ENDDO
178         ENDDO         ENDDO
       ENDDO  
       CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)  
 #endif /* ifndef SEAICE_EXTERNAL_FLUXES */  
179    
180  #ifdef SEAICE_TEST_ICE_STRESS_1        ELSE
181    C     else: useHB87StressCoupling=F
182    
183  C--   Compute ice-affected wind stress (interpolate to U/V-points)  C--   Compute ice-affected wind stress (interpolate to U/V-points)
184    C     by averaging wind stress and ice-ocean stress according to
185    C     ice cover
186        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
187         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
188          DO j=1,sNy          DO j=1,sNy
189           DO i=1,sNx           DO i=1,sNx
190            fuIce=HALF*( DWATN(I,J,bi,bj)+DWATN(I,J+1,bi,bj) )*(            fuIceLoc=HALF*( DWATN(I,J,bi,bj)+DWATN(I-1,J,bi,bj) )*
191       &         COSWAT *       &         COSWAT *
192       &         ( UICE(I,J,1,bi,bj)-GWATX(I,J,bi,bj) )       &         ( UICE(I,J,1,bi,bj)-uVel(I,J,1,bi,bj) )
193       &         - SINWAT* 0.5 _d 0 * (       &         - SIGN(SINWAT, _fCori(I,J,bi,bj)) * 0.5 _d 0 *
194       &          0.5 _d 0*(vIce(I  ,J  ,1,bi,bj)-GWATY(I  ,J  ,bi,bj)       &         ( DWATN(I  ,J,bi,bj) *
195       &                   +vIce(I-1,J  ,1,bi,bj)-GWATY(I-1,J  ,bi,bj))       &         0.5 _d 0*(vIce(I  ,J  ,1,bi,bj)-vVel(I  ,J  ,1,bi,bj)
196       &         +0.5 _d 0*(vIce(I  ,J+1,1,bi,bj)-GWATY(I  ,J+1,bi,bj)       &                  +vIce(I  ,J+1,1,bi,bj)-vVel(I  ,J+1,1,bi,bj))
197       &                   +vIce(I-1,J+1,1,bi,bj)-GWATY(I-1,J+1,bi,bj)) )       &         + DWATN(I-1,J,bi,bj) *
198         &         0.5 _d 0*(vIce(I-1,J  ,1,bi,bj)-vVel(I-1,J  ,1,bi,bj)
199         &                  +vIce(I-1,J+1,1,bi,bj)-vVel(I-1,J+1,1,bi,bj))
200       &         )       &         )
201            fvIce=HALF*( DWATN(I,J,bi,bj)+DWATN(I+1,J,bi,bj) )*(            fvIceLoc=HALF*( DWATN(I,J,bi,bj)+DWATN(I,J-1,bi,bj) )*
202       &         SINWAT *       &         COSWAT *
203       &         ( UICE(I,J,1,bi,bj)-GWATX(I,J,bi,bj) )       &         ( VICE(I,J,1,bi,bj)-vVel(I,J,1,bi,bj) )
204       &         + COSWAT * 0.5 _d 0 * (       &         + SIGN(SINWAT,  _fCori(I,J,bi,bj)) * 0.5 _d 0 *
205       &          0.5 _d 0*(uIce(I  ,J  ,1,bi,bj)-GWATY(I  ,J  ,bi,bj)       &         ( DWATN(I,J  ,bi,bj) *
206       &                   +uIce(I+1,J  ,1,bi,bj)-GWATX(I+1,J  ,bi,bj))       &         0.5 _d 0*(uIce(I  ,J  ,1,bi,bj)-uVel(I  ,J  ,1,bi,bj)
207       &         +0.5 _d 0*(uIce(I  ,J-1,1,bi,bj)-GWATY(I  ,J-1,bi,bj)       &                  +uIce(I+1,J  ,1,bi,bj)-uVel(I+1,J  ,1,bi,bj))
208       &                   +uIce(I+1,J-1,1,bi,bj)-GWATX(I+1,J-1,bi,bj)) )       &         + DWATN(I,J-1,bi,bj) *
209         &         0.5 _d 0*(uIce(I  ,J-1,1,bi,bj)-uVel(I  ,J-1,1,bi,bj)
210         &                  +uIce(I+1,J-1,1,bi,bj)-uVel(I+1,J-1,1,bi,bj))
211       &         )       &         )
212            fu(I,J,bi,bj)=(ONE-AREA(I,J,1,bi,bj))*fu(I,J,bi,bj)+            areaW = 0.5 _d 0 * (AREA(I,J,1,bi,bj) + AREA(I-1,J,1,bi,bj))
213       &         AREA(I,J,1,bi,bj)*fuIce       &         * SEAICEstressFactor
214            fv(I,J,bi,bj)=(ONE-AREA(I,J,1,bi,bj))*fv(I,J,bi,bj)+            areaS = 0.5 _d 0 * (AREA(I,J,1,bi,bj) + AREA(I,J-1,1,bi,bj))
215       &         AREA(I,J,1,bi,bj)*fvIce       &         * SEAICEstressFactor
216              fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)+areaW*fuIceLoc
217              fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)+areaS*fvIceLoc
218           ENDDO           ENDDO
219          ENDDO          ENDDO
220         ENDDO         ENDDO
221        ENDDO        ENDDO
222          ENDIF
223        CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)        CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)
224  #endif /* SEAICE_TEST_ICE_STRESS_1 */  
225  #endif /* not SEAICE_CGRID */  #endif /* SEAICE_CGRID */
226    
227        RETURN        RETURN
228        END        END
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.21
  ViewVC Help
Powered by ViewVC 1.1.22