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

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

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


Revision 1.2 - (show annotations) (download)
Mon Mar 6 13:27:40 2006 UTC (18 years, 2 months ago) by mlosch
Branch: MAIN
Changes since 1.1: +3 -3 lines
small bug fix

1 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_ocean_stress.F,v 1.1 2006/03/06 13:17:37 mlosch Exp $
2 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 #include "FFIELDS.h"
22 #include "SEAICE.h"
23 #include "SEAICE_PARAMS.h"
24
25 C === Routine arguments ===
26 C myTime - Simulation time
27 C myIter - Simulation timestep number
28 C myThid - Thread no. that called this routine.
29 _RL myTime
30 INTEGER myIter
31 INTEGER myThid
32 CML _RL COR_ICE (1-OLx:sNx+OLx,1-OLy:sNy+OLy, nSx,nSy)
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 _RL SINWIN, COSWIN, SINWAT, COSWAT
41 #ifdef SEAICE_TEST_ICE_STRESS_1
42 _RL fuIce, fvIce
43 #endif
44
45 c introduce turning angle (default is zero)
46 SINWIN=SIN(SEAICE_airTurnAngle*deg2rad)
47 COSWIN=COS(SEAICE_airTurnAngle*deg2rad)
48 SINWAT=SIN(SEAICE_waterTurnAngle*deg2rad)
49 COSWAT=COS(SEAICE_waterTurnAngle*deg2rad)
50
51 CML#ifdef SEAICE_ORIGINAL_BAD_ICE_STRESS
52 CMLC-- Following formulation is problematic and is no longer used.
53 CML#ifdef SEAICE_ALLOW_DYNAMICS
54 CML IF ( SEAICEuseDYNAMICS ) THEN
55 CMLC-- Compute ice-affected wind stress
56 CML DO bj=myByLo(myThid),myByHi(myThid)
57 CML DO bi=myBxLo(myThid),myBxHi(myThid)
58 CML DO j=1,sNy
59 CML DO i=1,sNx
60 CML WINDX(I,J,bi,bj)=DWATN(I,J,bi,bj)
61 CML & *(COSWAT*(GWATX(I,J,bi,bj)-UICE(I,J,1,bi,bj))
62 CML & -SINWAT*(GWATY(I,J,bi,bj)-VICEC(I,J,bi,bj)))
63 CML WINDY(I,J,bi,bj)=DWATN(I,J,bi,bj)
64 CML & *(SINWAT*(GWATX(I,J,bi,bj)-UICEC(I,J,bi,bj))
65 CML & +COSWAT*(GWATY(I,J,bi,bj)-VICE(I,J,1,bi,bj)))
66 CML WINDX(I,J,bi,bj)=WINDX(I,J,bi,bj)-( COR_ICE(I,J,bi,bj)
67 CML & *GWATY(I,J,bi,bj)-COR_ICE(I,J,bi,bj)*VICEC(I,J,bi,bj))
68 CML WINDY(I,J,bi,bj)=WINDY(I,J,bi,bj)-(-COR_ICE(I,J,bi,bj)
69 CML & *GWATX(I,J,bi,bj)+COR_ICE(I,J,bi,bj)*UICEC(I,J,bi,bj))
70 CML WINDX(I,J,bi,bj)=WINDX(I,J,bi,bj)-(UICE(I,J,1,bi,bj)
71 CML & -UICE(I,J,3,bi,bj))*AMASS(I,J,bi,bj)/SEAICE_DT*TWO
72 CML WINDY(I,J,bi,bj)=WINDY(I,J,bi,bj)-(VICE(I,J,1,bi,bj)
73 CML & -VICE(I,J,3,bi,bj))*AMASS(I,J,bi,bj)/SEAICE_DT*TWO
74 CML ENDDO
75 CML ENDDO
76 CML ENDDO
77 CML ENDDO
78 CML DO bj=myByLo(myThid),myByHi(myThid)
79 CML DO bi=myBxLo(myThid),myBxHi(myThid)
80 CML DO j=1,sNy
81 CML DO i=1,sNx
82 CML WINDX(I,J,bi,bj)=-WINDX(I,J,bi,bj)
83 CML WINDY(I,J,bi,bj)=-WINDY(I,J,bi,bj)
84 CML ENDDO
85 CML ENDDO
86 CML ENDDO
87 CML ENDDO
88 CML ENDIF
89 CML#endif /* SEAICE_ALLOW_DYNAMICS */
90 CML#endif /* SEAICE_ORIGINAL_BAD_ICE_STRESS */
91
92 C-- Update overlap regions
93 CALL EXCH_UV_XY_RL(WINDX, WINDY, .TRUE., myThid)
94
95 #ifndef SEAICE_EXTERNAL_FLUXES
96 C-- Interpolate wind stress (N/m^2) from South-West B-grid
97 C to South-West C-grid for forcing ocean model.
98 DO bj=myByLo(myThid),myByHi(myThid)
99 DO bi=myBxLo(myThid),myBxHi(myThid)
100 DO j=1,sNy
101 DO i=1,sNx
102 fu(I,J,bi,bj)=WINDX(I,J,bi,bj)
103 fv(I,J,bi,bj)=WINDY(I,J,bi,bj)
104 ENDDO
105 ENDDO
106 ENDDO
107 ENDDO
108 CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)
109 #endif /* ifndef SEAICE_EXTERNAL_FLUXES */
110
111 #ifdef SEAICE_TEST_ICE_STRESS_1
112 C-- Compute ice-affected wind stress (interpolate to U/V-points)
113 DO bj=myByLo(myThid),myByHi(myThid)
114 DO bi=myBxLo(myThid),myBxHi(myThid)
115 DO j=1,sNy
116 DO i=1,sNx
117 fuIce=HALF*( DWATN(I,J,bi,bj)+DWATN(I,J+1,bi,bj) )*(
118 & COSWAT *
119 & ( UICE(I,J,1,bi,bj)-GWATX(I,J,bi,bj) )
120 & - SINWAT* 0.5 _d 0 * (
121 & 0.5 _d 0*(vIce(I ,J ,1,bi,bj)-GWATY(I ,J ,bi,bj)
122 & +vIce(I-1,J ,1,bi,bj)-GWATY(I-1,J ,bi,bj))
123 & +0.5 _d 0*(vIce(I ,J+1,1,bi,bj)-GWATY(I ,J+1,bi,bj)
124 & +vIce(I-1,J+1,1,bi,bj)-GWATY(I-1,J+1,bi,bj)) )
125 & )
126 fvIce=HALF*( DWATN(I,J,bi,bj)+DWATN(I+1,J,bi,bj) )*(
127 & SINWAT *
128 & ( UICE(I,J,1,bi,bj)-GWATX(I,J,bi,bj) )
129 & + COSWAT * 0.5 _d 0 * (
130 & 0.5 _d 0*(uIce(I ,J ,1,bi,bj)-GWATX(I ,J ,bi,bj)
131 & +uIce(I+1,J ,1,bi,bj)-GWATX(I+1,J ,bi,bj))
132 & +0.5 _d 0*(uIce(I ,J-1,1,bi,bj)-GWATX(I ,J-1,bi,bj)
133 & +uIce(I+1,J-1,1,bi,bj)-GWATX(I+1,J-1,bi,bj)) )
134 & )
135 fu(I,J,bi,bj)=(ONE-AREA(I,J,1,bi,bj))*fu(I,J,bi,bj)+
136 & AREA(I,J,1,bi,bj)*fuIce
137 fv(I,J,bi,bj)=(ONE-AREA(I,J,1,bi,bj))*fv(I,J,bi,bj)+
138 & AREA(I,J,1,bi,bj)*fvIce
139 ENDDO
140 ENDDO
141 ENDDO
142 ENDDO
143 CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)
144 #endif /* SEAICE_TEST_ICE_STRESS_1 */
145 #endif /* not SEAICE_CGRID */
146
147 RETURN
148 END

  ViewVC Help
Powered by ViewVC 1.1.22