/[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.18 - (hide annotations) (download)
Tue May 15 14:32:56 2007 UTC (17 years ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint59e, checkpoint59d, checkpoint59g, checkpoint59f, checkpoint59c, checkpoint59b, checkpoint59i, checkpoint59h, checkpoint59j
Changes since 1.17: +10 -5 lines 
  - get rid of three more exchanges seaice_calc_viscosities (required
    different loop boundaries in seaice_lsr and seaice_ocean_stress)
  - bug fix in seaice_ocean_stress (update lab_sea)
  - reduce loop boundaries to necessary minimum also in seaice_evp
  - go back to the old way of computing deltaC/Z as the adjoint likes that
    better
1 mlosch 1.18 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_ocean_stress.F,v 1.17 2007/05/14 14:29:49 jmc 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 mlosch 1.5 IF ( useHB87StressCoupling ) THEN
60     C
61     C use an intergral over ice and ocean surface layer to define
62     C surface stresses on ocean following Hibler and Bryan (1987, JPO)
63     C
64 mlosch 1.14 C recompute strain rates, viscosities, etc. from updated ice velocities
65     IF ( .NOT. SEAICEuseEVP ) THEN
66 mlosch 1.16 C only for EVP we already have the stress components otherwise we need
67     C to recompute them here
68 mlosch 1.14 CALL SEAICE_CALC_STRAINRATES(
69     I uIce(1-Olx,1-Oly,1,1,1), vIce(1-Olx,1-Oly,1,1,1),
70     O e11, e22, e12,
71     I myThid )
72    
73     CALL SEAICE_CALC_VISCOSITIES(
74     I e11, e22, e12, zMin, zMax, hEffM, press0,
75     O eta, zeta, press,
76     I myThid )
77     ENDIF
78 mlosch 1.5 C re-compute internal stresses with updated ice velocities
79     DO bj=myByLo(myThid),myByHi(myThid)
80     DO bi=myBxLo(myThid),myBxHi(myThid)
81 mlosch 1.14 IF ( .NOT. SEAICEuseEVP ) THEN
82     C only for EVP we already have computed the stress divergences, for
83     C anything else we have to do it here
84     DO j=1-Oly,sNy+Oly
85     DO i=1-Olx,sNx+Olx
86     sig11(I,J) = 0. _d 0
87     sig22(I,J) = 0. _d 0
88     sig12(I,J) = 0. _d 0
89     ENDDO
90 mlosch 1.5 ENDDO
91 mlosch 1.14
92 mlosch 1.18 DO j=0,sNy
93     DO i=0,sNx
94 mlosch 1.14 eplus = e11(I,J,bi,bj) + e22(I,J,bi,bj)
95     eminus= e11(I,J,bi,bj) - e22(I,J,bi,bj)
96     sig11(I,J) = zeta(I,J,bi,bj)*eplus + eta(I,J,bi,bj)*eminus
97     & - 0.5 _d 0 * PRESS(I,J,bi,bj)
98     sig22(I,J) = zeta(I,J,bi,bj)*eplus - eta(I,J,bi,bj)*eminus
99     & - 0.5 _d 0 * PRESS(I,J,bi,bj)
100 mlosch 1.18 ENDDO
101     ENDDO
102    
103     DO j=1,sNy+1
104     DO i=1,sNx+1
105 mlosch 1.14 sig12(I,J) = 2. _d 0 * e12(I,J,bi,bj) *
106     & ( 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) )
108     & /MAX(1. _d 0,
109     & 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))
111     ENDDO
112     ENDDO
113     C evaluate divergence of stress and apply to forcing
114     DO J=1,sNy
115     DO I=1,sNx
116     FX = ( sig11(I ,J ) * _dyF(I ,J ,bi,bj)
117     & - sig11(I-1,J ) * _dyF(I-1,J ,bi,bj)
118     & + sig12(I ,J+1) * _dxV(I ,J+1,bi,bj)
119     & - sig12(I ,J ) * _dxV(I ,J ,bi,bj)
120     & ) * recip_rAw(I,J,bi,bj)
121     & -
122     & ( sig12(I,J) + sig12(I,J+1) )
123     & * _tanPhiAtU(I,J,bi,bj) * recip_rSphere
124     & +
125     & ( sig22(I,J) + sig22(I-1,J) ) * 0.5 _d 0
126     & * _tanPhiAtU(I,J,bi,bj) * recip_rSphere
127     C one metric term missing for general curvilinear coordinates
128     FY = ( sig22(I ,J ) * _dxF(I ,J ,bi,bj)
129     & - sig22(I ,J-1) * _dxF(I ,J-1,bi,bj)
130     & + sig12(I+1,J ) * _dyU(I+1,J ,bi,bj)
131     & - sig12(I ,J ) * _dyU(I ,J ,bi,bj)
132     & ) * recip_rAs(I,J,bi,bj)
133     & -
134     & ( sig22(I,J) + sig22(I,J-1) ) * 0.5 _d 0
135     & * _tanPhiAtV(I,J,bi,bj) * recip_rSphere
136     C two metric terms missing for general curvilinear coordinates
137     C average wind stress over ice and ocean and apply averaged wind
138     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))
140     & * SEAICEstressFactor
141     areaS = 0.5 * (AREA(I,J,1,bi,bj) + AREA(I,J-1,1,bi,bj))
142     & * SEAICEstressFactor
143     fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)
144     & + areaW*taux(I,J,bi,bj)
145     & + FX * SEAICEstressFactor
146     fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)
147     & + areaS*tauy(I,J,bi,bj)
148     & + FY * SEAICEstressFactor
149     C save stress divergence for later
150 mlosch 1.16 #ifdef SEAICE_ALLOW_EVP
151 mlosch 1.14 stressDivergenceX(I,J,bi,bj) = FX
152     stressDivergenceY(I,J,bi,bj) = FY
153 mlosch 1.16 #endif /* SEAICE_ALLOW_EVP */
154 mlosch 1.14 ENDDO
155     ENDDO
156     ELSE
157 mlosch 1.16 #ifdef SEAICE_ALLOW_EVP
158 mlosch 1.14 DO J=1,sNy
159     DO I=1,sNx
160 mlosch 1.5 C average wind stress over ice and ocean and apply averaged wind
161     C stress and internal ice stresses to surface layer of ocean
162 mlosch 1.14 areaW = 0.5 * (AREA(I,J,1,bi,bj) + AREA(I-1,J,1,bi,bj))
163     & * SEAICEstressFactor
164     areaS = 0.5 * (AREA(I,J,1,bi,bj) + AREA(I,J-1,1,bi,bj))
165     & * SEAICEstressFactor
166     fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)
167     & + areaW*taux(I,J,bi,bj)
168     & + stressDivergenceX(I,J,bi,bj) * SEAICEstressFactor
169     fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)
170     & + areaS*tauy(I,J,bi,bj)
171     & + stressDivergenceY(I,J,bi,bj) * SEAICEstressFactor
172     ENDDO
173     ENDDO
174 mlosch 1.16 #endif /* SEAICE_ALLOW_EVP */
175 mlosch 1.14 ENDIF
176 mlosch 1.5 ENDDO
177     ENDDO
178 jmc 1.17
179 mlosch 1.5 ELSE
180 jmc 1.17 C else: useHB87StressCoupling=F
181 mlosch 1.5
182     C-- Compute ice-affected wind stress (interpolate to U/V-points)
183     C by averaging wind stress and ice-ocean stress according to
184     C ice cover
185 mlosch 1.1 DO bj=myByLo(myThid),myByHi(myThid)
186     DO bi=myBxLo(myThid),myBxHi(myThid)
187     DO j=1,sNy
188     DO i=1,sNx
189 mlosch 1.18 fuIceLoc=HALF*( DWATN(I,J,bi,bj)+DWATN(I-1,J,bi,bj) )*
190 mlosch 1.1 & COSWAT *
191     & ( UICE(I,J,1,bi,bj)-GWATX(I,J,bi,bj) )
192 mlosch 1.6 & - SIGN(SINWAT, _fCori(I,J,bi,bj)) * 0.5 _d 0 *
193     & ( DWATN(I ,J,bi,bj) *
194     & 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))
196     & + DWATN(I-1,J,bi,bj) *
197     & 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))
199 mlosch 1.1 & )
200 mlosch 1.18 fvIceLoc=HALF*( DWATN(I,J,bi,bj)+DWATN(I,J-1,bi,bj) )*
201 mlosch 1.6 & COSWAT *
202     & ( VICE(I,J,1,bi,bj)-GWATY(I,J,bi,bj) )
203     & + SIGN(SINWAT, _fCori(I,J,bi,bj)) * 0.5 _d 0 *
204     & ( DWATN(I,J ,bi,bj) *
205     & 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))
207     & + DWATN(I,J-1,bi,bj) *
208     & 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))
210 mlosch 1.1 & )
211 mlosch 1.4 areaW = 0.5 _d 0 * (AREA(I,J,1,bi,bj) + AREA(I-1,J,1,bi,bj))
212 mlosch 1.9 & * SEAICEstressFactor
213 mlosch 1.4 areaS = 0.5 _d 0 * (AREA(I,J,1,bi,bj) + AREA(I,J-1,1,bi,bj))
214 mlosch 1.9 & * SEAICEstressFactor
215 mlosch 1.11 fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)+areaW*fuIceLoc
216     fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)+areaS*fvIceLoc
217 mlosch 1.1 ENDDO
218     ENDDO
219     ENDDO
220     ENDDO
221 mlosch 1.5 ENDIF
222 mlosch 1.1 CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)
223 mlosch 1.3
224 jmc 1.17 #endif /* SEAICE_CGRID */
225 mlosch 1.1
226     RETURN
227     END
  ViewVC Help
Powered by ViewVC 1.1.22