/[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.4 by mlosch, Thu Mar 9 20:22:40 2006 UTC revision 1.23 by mlosch, Fri May 29 10:18:03 2009 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  SINWAT, COSWAT        _RL  SINWAT, COSWAT, SINWIN, COSWIN
42        _RL  fuIce, fvIce        _RL  fuIceLoc, fvIceLoc, FX, FY
43        _RL  areaW, areaS        _RL  areaW, areaS
44    
45          _RL sig11       (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
46          _RL sig22       (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
47          _RL sig12       (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
48          _RL eplus, eminus
49    
50  c     introduce turning angle (default is zero)  c     introduce turning angle (default is zero)
51        SINWAT=SIN(SEAICE_waterTurnAngle*deg2rad)        SINWAT=SIN(SEAICE_waterTurnAngle*deg2rad)
52        COSWAT=COS(SEAICE_waterTurnAngle*deg2rad)        COSWAT=COS(SEAICE_waterTurnAngle*deg2rad)
53          SINWIN=SIN(SEAICE_airTurnAngle*deg2rad)
54          COSWIN=COS(SEAICE_airTurnAngle*deg2rad)
55    
56  C--   Update overlap regions        IF ( useHB87StressCoupling ) THEN
57        CALL EXCH_UV_XY_RL(WINDX, WINDY, .TRUE., myThid)  C
58    C     use an intergral over ice and ocean surface layer to define
59  #ifndef SEAICE_EXTERNAL_FLUXES  C     surface stresses on ocean following Hibler and Bryan (1987, JPO)
60  C--   Interpolate wind stress (N/m^2) from C-points of C-grid  C
61  C     to U and V points of C-grid for forcing the ocean model.  C     recompute strain rates, viscosities, etc. from updated ice velocities
62        DO bj=myByLo(myThid),myByHi(myThid)         IF ( .NOT. SEAICEuseEVP ) THEN
63         DO bi=myBxLo(myThid),myBxHi(myThid)  C     only for EVP we already have the stress components otherwise we need
64          DO j=1,sNy  C     to recompute them here
65           DO i=1,sNx          CALL SEAICE_CALC_STRAINRATES(
66            fu(I,J,bi,bj)=0.5*(WINDX(I,J,bi,bj) + WINDX(I-1,J,bi,bj))       I       uIce, vIce,
67            fv(I,J,bi,bj)=0.5*(WINDY(I,J,bi,bj) + WINDY(I,J-1,bi,bj))       O       e11, e22, e12,
68           ENDDO       I       3, 3, myTime, myIter, myThid )
69    
70            CALL SEAICE_CALC_VISCOSITIES(
71         I       e11, e22, e12, zMin, zMax, hEffM, press0,
72         O       eta, zeta, press,
73         I       3, myTime, myIter, myThid )
74           ENDIF
75    C     re-compute internal stresses with updated ice velocities
76           DO bj=myByLo(myThid),myByHi(myThid)
77            DO bi=myBxLo(myThid),myBxHi(myThid)
78             IF ( .NOT. SEAICEuseEVP ) THEN
79    C     only for EVP we already have computed the stress divergences, for
80    C     anything else we have to do it here
81              DO j=1-Oly,sNy+Oly
82               DO i=1-Olx,sNx+Olx
83                sig11(I,J) = 0. _d 0
84                sig22(I,J) = 0. _d 0
85                sig12(I,J) = 0. _d 0
86               ENDDO
87              ENDDO
88    
89              DO j=0,sNy
90               DO i=0,sNx
91                eplus = e11(I,J,bi,bj) + e22(I,J,bi,bj)
92                eminus= e11(I,J,bi,bj) - e22(I,J,bi,bj)
93                sig11(I,J) = zeta(I,J,bi,bj)*eplus + eta(I,J,bi,bj)*eminus
94         &           - 0.5 _d 0 * PRESS(I,J,bi,bj)
95                sig22(I,J) = zeta(I,J,bi,bj)*eplus - eta(I,J,bi,bj)*eminus
96         &           - 0.5 _d 0 * PRESS(I,J,bi,bj)
97               ENDDO
98              ENDDO
99    
100              DO j=1,sNy+1
101               DO i=1,sNx+1
102                sig12(I,J) = 2. _d 0 * e12(I,J,bi,bj) *
103         &           ( eta(I,J  ,bi,bj) + eta(I-1,J  ,bi,bj)
104         &           + eta(I,J-1,bi,bj) + eta(I-1,J-1,bi,bj) )
105         &           /MAX(1. _d 0,
106         &             hEffM(I,J  ,bi,bj) + hEffM(I-1,J  ,bi,bj)
107         &           + hEffM(I,J-1,bi,bj) + hEffM(I-1,J-1,bi,bj))
108               ENDDO
109              ENDDO
110    C     evaluate divergence of stress and apply to forcing
111              DO J=1,sNy
112               DO I=1,sNx
113                FX = ( sig11(I  ,J  ) * _dyF(I  ,J  ,bi,bj)
114         &           - sig11(I-1,J  ) * _dyF(I-1,J  ,bi,bj)
115         &           + sig12(I  ,J+1) * _dxV(I  ,J+1,bi,bj)
116         &           - sig12(I  ,J  ) * _dxV(I  ,J  ,bi,bj)
117         &           ) * recip_rAw(I,J,bi,bj)
118                FY = ( sig22(I  ,J  ) * _dxF(I  ,J  ,bi,bj)
119         &           - sig22(I  ,J-1) * _dxF(I  ,J-1,bi,bj)
120         &           + sig12(I+1,J  ) * _dyU(I+1,J  ,bi,bj)
121         &           - sig12(I  ,J  ) * _dyU(I  ,J  ,bi,bj)
122         &           ) * recip_rAs(I,J,bi,bj)
123    C     average wind stress over ice and ocean and apply averaged wind
124    C     stress and internal ice stresses to surface layer of ocean
125                areaW = 0.5 * (AREA(I,J,1,bi,bj) + AREA(I-1,J,1,bi,bj))
126         &           * SEAICEstressFactor
127                areaS = 0.5 * (AREA(I,J,1,bi,bj) + AREA(I,J-1,1,bi,bj))
128         &           * SEAICEstressFactor
129                fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)
130         &           + areaW*taux(I,J,bi,bj)
131         &           + FX * SEAICEstressFactor
132                fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)
133         &           + areaS*tauy(I,J,bi,bj)
134         &           + FY * SEAICEstressFactor
135    C     save stress divergence for later
136    #ifdef SEAICE_ALLOW_EVP
137                stressDivergenceX(I,J,bi,bj) = FX
138                stressDivergenceY(I,J,bi,bj) = FY
139    #endif /* SEAICE_ALLOW_EVP */
140               ENDDO
141              ENDDO
142             ELSE
143    #ifdef SEAICE_ALLOW_EVP
144              DO J=1,sNy
145               DO I=1,sNx
146    C     average wind stress over ice and ocean and apply averaged wind
147    C     stress and internal ice stresses to surface layer of ocean
148                areaW = 0.5 * (AREA(I,J,1,bi,bj) + AREA(I-1,J,1,bi,bj))
149         &           * SEAICEstressFactor
150                areaS = 0.5 * (AREA(I,J,1,bi,bj) + AREA(I,J-1,1,bi,bj))
151         &           * SEAICEstressFactor
152                fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)
153         &           + areaW*taux(I,J,bi,bj)
154         &           + stressDivergenceX(I,J,bi,bj) * SEAICEstressFactor
155                fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)
156         &           + areaS*tauy(I,J,bi,bj)
157         &           + stressDivergenceY(I,J,bi,bj) * SEAICEstressFactor
158               ENDDO
159              ENDDO
160    #endif /* SEAICE_ALLOW_EVP */
161             ENDIF
162          ENDDO          ENDDO
163         ENDDO         ENDDO
164        ENDDO  
165  #endif /* ifndef SEAICE_EXTERNAL_FLUXES */        ELSE
166    C     else: useHB87StressCoupling=F
167    
168  C--   Compute ice-affected wind stress (interpolate to U/V-points)  C--   Compute ice-affected wind stress (interpolate to U/V-points)
169    C     by averaging wind stress and ice-ocean stress according to
170    C     ice cover
171        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
172         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
173          DO j=1,sNy          DO j=1,sNy
174           DO i=1,sNx           DO i=1,sNx
175            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) )*
176       &         COSWAT *       &         COSWAT *
177       &         ( UICE(I,J,1,bi,bj)-GWATX(I,J,bi,bj) )       &         ( UICE(I,J,1,bi,bj)-uVel(I,J,1,bi,bj) )
178       &         - SINWAT* 0.5 _d 0 * (       &         - SIGN(SINWAT, _fCori(I,J,bi,bj)) * 0.5 _d 0 *
179       &          0.5 _d 0*(vIce(I  ,J  ,1,bi,bj)-GWATY(I  ,J  ,bi,bj)       &         ( DWATN(I  ,J,bi,bj) *
180       &                   +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)
181       &         +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))
182       &                   +vIce(I-1,J+1,1,bi,bj)-GWATY(I-1,J+1,bi,bj)) )       &         + DWATN(I-1,J,bi,bj) *
183         &         0.5 _d 0*(vIce(I-1,J  ,1,bi,bj)-vVel(I-1,J  ,1,bi,bj)
184         &                  +vIce(I-1,J+1,1,bi,bj)-vVel(I-1,J+1,1,bi,bj))
185       &         )       &         )
186            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) )*
187       &         SINWAT *       &         COSWAT *
188       &         ( UICE(I,J,1,bi,bj)-GWATX(I,J,bi,bj) )       &         ( VICE(I,J,1,bi,bj)-vVel(I,J,1,bi,bj) )
189       &         + COSWAT * 0.5 _d 0 * (       &         + SIGN(SINWAT,  _fCori(I,J,bi,bj)) * 0.5 _d 0 *
190       &          0.5 _d 0*(uIce(I  ,J  ,1,bi,bj)-GWATX(I  ,J  ,bi,bj)       &         ( DWATN(I,J  ,bi,bj) *
191       &                   +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)
192       &         +0.5 _d 0*(uIce(I  ,J-1,1,bi,bj)-GWATX(I  ,J-1,bi,bj)       &                  +uIce(I+1,J  ,1,bi,bj)-uVel(I+1,J  ,1,bi,bj))
193       &                   +uIce(I+1,J-1,1,bi,bj)-GWATX(I+1,J-1,bi,bj)) )       &         + DWATN(I,J-1,bi,bj) *
194         &         0.5 _d 0*(uIce(I  ,J-1,1,bi,bj)-uVel(I  ,J-1,1,bi,bj)
195         &                  +uIce(I+1,J-1,1,bi,bj)-uVel(I+1,J-1,1,bi,bj))
196       &         )       &         )
197            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))
198         &         * SEAICEstressFactor
199            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))
200            fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)+areaW*fuIce       &         * SEAICEstressFactor
201            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
202              fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)+areaS*fvIceLoc
203           ENDDO           ENDDO
204          ENDDO          ENDDO
205         ENDDO         ENDDO
206        ENDDO        ENDDO
207          ENDIF
208        CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)        CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)
209    
210  #endif /* not SEAICE_CGRID */  #endif /* SEAICE_CGRID */
211    
212        RETURN        RETURN
213        END        END
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.23
  ViewVC Help
Powered by ViewVC 1.1.22