/[MITgcm]/MITgcm/pkg/seaice/seaice_ocean_stress.F
ViewVC logotype

Annotation of /MITgcm/pkg/seaice/seaice_ocean_stress.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph

Revision 1.14 - (hide annotations) (download)
Tue Apr 24 11:23:26 2007 UTC (17 years, 1 month ago) by mlosch
Branch: MAIN
Changes since 1.13: +103 -159 lines 
  - modify Hibler+Bryan (1987) stress coupling so that it reuses information
    of evp, complete rewrite for lsr case (latter is does not yet work
    in all cases)
1 mlosch 1.14 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_ocean_stress.F,v 1.13 2007/04/20 18:29:58 mlosch Exp $
2 mlosch 1.1 C $Name: $
3    
4     #include "SEAICE_OPTIONS.h"
5    
6     CStartOfInterface
7     SUBROUTINE SEAICE_OCEAN_STRESS(
8     I myTime, myIter, myThid )
9     C /==========================================================\
10     C | SUBROUTINE SEAICE_OCEAN_STRESS |
11     C | o Calculate ocean surface stresses |
12     C | - C-grid version |
13     C |==========================================================|
14     C \==========================================================/
15     IMPLICIT NONE
16    
17     C === Global variables ===
18     #include "SIZE.h"
19     #include "EEPARAMS.h"
20     #include "PARAMS.h"
21 mlosch 1.5 #include "GRID.h"
22 mlosch 1.1 #include "FFIELDS.h"
23     #include "SEAICE.h"
24     #include "SEAICE_PARAMS.h"
25    
26     C === Routine arguments ===
27     C myTime - Simulation time
28     C myIter - Simulation timestep number
29     C myThid - Thread no. that called this routine.
30     _RL myTime
31     INTEGER myIter
32     INTEGER myThid
33     CEndOfInterface
34    
35     #ifdef SEAICE_CGRID
36     C === Local variables ===
37     C i,j,bi,bj - Loop counters
38    
39     INTEGER i, j, bi, bj
40 mlosch 1.5 _RL SINWAT, COSWAT, SINWIN, COSWIN
41 mlosch 1.11 _RL fuIceLoc, fvIceLoc, FX, FY
42 mlosch 1.4 _RL areaW, areaS
43 mlosch 1.1
44 mlosch 1.13 _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 mlosch 1.5 _RL press (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
48 mlosch 1.14 _RL sig11 (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
49     _RL sig22 (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
50     _RL sig12 (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
51     _RL eplus, eminus
52 mlosch 1.5
53 mlosch 1.1 c introduce turning angle (default is zero)
54     SINWAT=SIN(SEAICE_waterTurnAngle*deg2rad)
55     COSWAT=COS(SEAICE_waterTurnAngle*deg2rad)
56 mlosch 1.5 SINWIN=SIN(SEAICE_airTurnAngle*deg2rad)
57     COSWIN=COS(SEAICE_airTurnAngle*deg2rad)
58 mlosch 1.1
59     C-- Update overlap regions
60     CALL EXCH_UV_XY_RL(WINDX, WINDY, .TRUE., myThid)
61    
62     #ifndef SEAICE_EXTERNAL_FLUXES
63 mlosch 1.3 C-- Interpolate wind stress (N/m^2) from C-points of C-grid
64     C to U and V points of C-grid for forcing the ocean model.
65 mlosch 1.1 DO bj=myByLo(myThid),myByHi(myThid)
66     DO bi=myBxLo(myThid),myBxHi(myThid)
67     DO j=1,sNy
68     DO i=1,sNx
69 mlosch 1.3 fu(I,J,bi,bj)=0.5*(WINDX(I,J,bi,bj) + WINDX(I-1,J,bi,bj))
70     fv(I,J,bi,bj)=0.5*(WINDY(I,J,bi,bj) + WINDY(I,J-1,bi,bj))
71 mlosch 1.1 ENDDO
72     ENDDO
73     ENDDO
74     ENDDO
75     #endif /* ifndef SEAICE_EXTERNAL_FLUXES */
76    
77 mlosch 1.5 IF ( useHB87StressCoupling ) THEN
78     C
79     C use an intergral over ice and ocean surface layer to define
80     C surface stresses on ocean following Hibler and Bryan (1987, JPO)
81     C
82 mlosch 1.14 C recompute strain rates, viscosities, etc. from updated ice velocities
83     IF ( .NOT. SEAICEuseEVP ) THEN
84     C we already have the stress components and do not need to recompute them
85     CALL SEAICE_CALC_STRAINRATES(
86     I uIce(1-Olx,1-Oly,1,1,1), vIce(1-Olx,1-Oly,1,1,1),
87     O e11, e22, e12,
88     I myThid )
89    
90     CALL SEAICE_CALC_VISCOSITIES(
91     I e11, e22, e12, zMin, zMax, hEffM, press0,
92     O eta, zeta, press,
93     I myThid )
94     ENDIF
95 mlosch 1.5 C re-compute internal stresses with updated ice velocities
96     DO bj=myByLo(myThid),myByHi(myThid)
97     DO bi=myBxLo(myThid),myBxHi(myThid)
98 mlosch 1.14 IF ( .NOT. SEAICEuseEVP ) THEN
99     C only for EVP we already have computed the stress divergences, for
100     C anything else we have to do it here
101     DO j=1-Oly,sNy+Oly
102     DO i=1-Olx,sNx+Olx
103     sig11(I,J) = 0. _d 0
104     sig22(I,J) = 0. _d 0
105     sig12(I,J) = 0. _d 0
106     ENDDO
107 mlosch 1.5 ENDDO
108 mlosch 1.14
109     DO j=1-Oly+1,sNy+Oly-1
110     DO i=1-Olx+1,sNx+Olx-1
111     eplus = e11(I,J,bi,bj) + e22(I,J,bi,bj)
112     eminus= e11(I,J,bi,bj) - e22(I,J,bi,bj)
113     sig11(I,J) = zeta(I,J,bi,bj)*eplus + eta(I,J,bi,bj)*eminus
114     & - 0.5 _d 0 * PRESS(I,J,bi,bj)
115     sig22(I,J) = zeta(I,J,bi,bj)*eplus - eta(I,J,bi,bj)*eminus
116     & - 0.5 _d 0 * PRESS(I,J,bi,bj)
117     sig12(I,J) = 2. _d 0 * e12(I,J,bi,bj) *
118     & ( eta(I,J ,bi,bj) + eta(I-1,J ,bi,bj)
119     & + eta(I,J-1,bi,bj) + eta(I-1,J-1,bi,bj) )
120     & /MAX(1. _d 0,
121     & hEffM(I,J ,bi,bj) + hEffM(I-1,J ,bi,bj)
122     & + hEffM(I,J-1,bi,bj) + hEffM(I-1,J-1,bi,bj))
123     ENDDO
124     ENDDO
125     C evaluate divergence of stress and apply to forcing
126     DO J=1,sNy
127     DO I=1,sNx
128     FX = ( sig11(I ,J ) * _dyF(I ,J ,bi,bj)
129     & - sig11(I-1,J ) * _dyF(I-1,J ,bi,bj)
130     & + sig12(I ,J+1) * _dxV(I ,J+1,bi,bj)
131     & - sig12(I ,J ) * _dxV(I ,J ,bi,bj)
132     & ) * recip_rAw(I,J,bi,bj)
133     & -
134     & ( sig12(I,J) + sig12(I,J+1) )
135     & * _tanPhiAtU(I,J,bi,bj) * recip_rSphere
136     & +
137     & ( sig22(I,J) + sig22(I-1,J) ) * 0.5 _d 0
138     & * _tanPhiAtU(I,J,bi,bj) * recip_rSphere
139     C one metric term missing for general curvilinear coordinates
140     FY = ( sig22(I ,J ) * _dxF(I ,J ,bi,bj)
141     & - sig22(I ,J-1) * _dxF(I ,J-1,bi,bj)
142     & + sig12(I+1,J ) * _dyU(I+1,J ,bi,bj)
143     & - sig12(I ,J ) * _dyU(I ,J ,bi,bj)
144     & ) * recip_rAs(I,J,bi,bj)
145     & -
146     & ( sig22(I,J) + sig22(I,J-1) ) * 0.5 _d 0
147     & * _tanPhiAtV(I,J,bi,bj) * recip_rSphere
148     C two metric terms missing for general curvilinear coordinates
149     C average wind stress over ice and ocean and apply averaged wind
150     C stress and internal ice stresses to surface layer of ocean
151     areaW = 0.5 * (AREA(I,J,1,bi,bj) + AREA(I-1,J,1,bi,bj))
152     & * SEAICEstressFactor
153     areaS = 0.5 * (AREA(I,J,1,bi,bj) + AREA(I,J-1,1,bi,bj))
154     & * SEAICEstressFactor
155     fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)
156     & + areaW*taux(I,J,bi,bj)
157     & + FX * SEAICEstressFactor
158     fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)
159     & + areaS*tauy(I,J,bi,bj)
160     & + FY * SEAICEstressFactor
161     C save stress divergence for later
162     #ifdef ALLOW_EVP
163     stressDivergenceX(I,J,bi,bj) = FX
164     stressDivergenceY(I,J,bi,bj) = FY
165     #endif
166     ENDDO
167     ENDDO
168     ELSE
169     DO J=1,sNy
170     DO I=1,sNx
171 mlosch 1.5 C average wind stress over ice and ocean and apply averaged wind
172     C stress and internal ice stresses to surface layer of ocean
173 mlosch 1.14 areaW = 0.5 * (AREA(I,J,1,bi,bj) + AREA(I-1,J,1,bi,bj))
174     & * SEAICEstressFactor
175     areaS = 0.5 * (AREA(I,J,1,bi,bj) + AREA(I,J-1,1,bi,bj))
176     & * SEAICEstressFactor
177     fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)
178     & + areaW*taux(I,J,bi,bj)
179     & + stressDivergenceX(I,J,bi,bj) * SEAICEstressFactor
180     fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)
181     & + areaS*tauy(I,J,bi,bj)
182     & + stressDivergenceY(I,J,bi,bj) * SEAICEstressFactor
183     ENDDO
184     ENDDO
185     ENDIF
186 mlosch 1.5 ENDDO
187     ENDDO
188     ELSE
189    
190     C-- Compute ice-affected wind stress (interpolate to U/V-points)
191     C by averaging wind stress and ice-ocean stress according to
192     C ice cover
193 mlosch 1.1 DO bj=myByLo(myThid),myByHi(myThid)
194     DO bi=myBxLo(myThid),myBxHi(myThid)
195     DO j=1,sNy
196     DO i=1,sNx
197 mlosch 1.11 fuIceLoc=HALF*( DWATN(I,J,bi,bj)+DWATN(I,J+1,bi,bj) )*
198 mlosch 1.1 & COSWAT *
199     & ( UICE(I,J,1,bi,bj)-GWATX(I,J,bi,bj) )
200 mlosch 1.6 & - SIGN(SINWAT, _fCori(I,J,bi,bj)) * 0.5 _d 0 *
201     & ( DWATN(I ,J,bi,bj) *
202     & 0.5 _d 0*(vIce(I ,J ,1,bi,bj)-GWATY(I ,J ,bi,bj)
203     & +vIce(I ,J+1,1,bi,bj)-GWATY(I ,J+1,bi,bj))
204     & + DWATN(I-1,J,bi,bj) *
205     & 0.5 _d 0*(vIce(I-1,J ,1,bi,bj)-GWATY(I-1,J ,bi,bj)
206     & +vIce(I-1,J+1,1,bi,bj)-GWATY(I-1,J+1,bi,bj))
207 mlosch 1.1 & )
208 mlosch 1.11 fvIceLoc=HALF*( DWATN(I,J,bi,bj)+DWATN(I+1,J,bi,bj) )*
209 mlosch 1.6 & COSWAT *
210     & ( VICE(I,J,1,bi,bj)-GWATY(I,J,bi,bj) )
211     & + SIGN(SINWAT, _fCori(I,J,bi,bj)) * 0.5 _d 0 *
212     & ( DWATN(I,J ,bi,bj) *
213     & 0.5 _d 0*(uIce(I ,J ,1,bi,bj)-GWATX(I ,J ,bi,bj)
214     & +uIce(I+1,J ,1,bi,bj)-GWATX(I+1,J ,bi,bj))
215     & + DWATN(I,J-1,bi,bj) *
216     & 0.5 _d 0*(uIce(I ,J-1,1,bi,bj)-GWATX(I ,J-1,bi,bj)
217     & +uIce(I+1,J-1,1,bi,bj)-GWATX(I+1,J-1,bi,bj))
218 mlosch 1.1 & )
219 mlosch 1.4 areaW = 0.5 _d 0 * (AREA(I,J,1,bi,bj) + AREA(I-1,J,1,bi,bj))
220 mlosch 1.9 & * SEAICEstressFactor
221 mlosch 1.4 areaS = 0.5 _d 0 * (AREA(I,J,1,bi,bj) + AREA(I,J-1,1,bi,bj))
222 mlosch 1.9 & * SEAICEstressFactor
223 mlosch 1.11 fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)+areaW*fuIceLoc
224     fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)+areaS*fvIceLoc
225 mlosch 1.1 ENDDO
226     ENDDO
227     ENDDO
228     ENDDO
229 mlosch 1.5 ENDIF
230 mlosch 1.1 CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)
231 mlosch 1.3
232 mlosch 1.1 #endif /* not SEAICE_CGRID */
233    
234     RETURN
235     END
  ViewVC Help
Powered by ViewVC 1.1.22