/[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.2 by mlosch, Mon Mar 6 13:27:40 2006 UTC revision 1.27 by jmc, Thu Sep 10 16:05:26 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    C     kSrf      - vertical index of surface layer
40        INTEGER i, j, bi, bj        INTEGER i, j, bi, bj
41        _RL  SINWIN, COSWIN, SINWAT, COSWAT        INTEGER kSrf
42  #ifdef SEAICE_TEST_ICE_STRESS_1        _RL  COSWAT
43        _RL  fuIce, fvIce        _RS  SINWAT
44  #endif        _RL  fuIceLoc, fvIceLoc
45          _RL  areaW, areaS
46  c     introduce turning angle (default is zero)  
47        SINWIN=SIN(SEAICE_airTurnAngle*deg2rad)  C     surrface level
48        COSWIN=COS(SEAICE_airTurnAngle*deg2rad)        kSrf = 1
49    C     introduce turning angle (default is zero)
50        SINWAT=SIN(SEAICE_waterTurnAngle*deg2rad)        SINWAT=SIN(SEAICE_waterTurnAngle*deg2rad)
51        COSWAT=COS(SEAICE_waterTurnAngle*deg2rad)        COSWAT=COS(SEAICE_waterTurnAngle*deg2rad)
52    
53  CML#ifdef SEAICE_ORIGINAL_BAD_ICE_STRESS        IF ( useHB87StressCoupling ) THEN
54  CMLC--   Following formulation is problematic and is no longer used.  C
55  CML#ifdef SEAICE_ALLOW_DYNAMICS  C     use an intergral over ice and ocean surface layer to define
56  CML      IF ( SEAICEuseDYNAMICS ) THEN  C     surface stresses on ocean following Hibler and Bryan (1987, JPO)
57  CMLC--   Compute ice-affected wind stress  C
58  CML       DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
59  CML        DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
60  CML         DO j=1,sNy           DO J=1,sNy
61  CML          DO i=1,sNx            DO I=1,sNx
62  CML           WINDX(I,J,bi,bj)=DWATN(I,J,bi,bj)  C     average wind stress over ice and ocean and apply averaged wind
63  CML     &          *(COSWAT*(GWATX(I,J,bi,bj)-UICE(I,J,1,bi,bj))  C     stress and internal ice stresses to surface layer of ocean
64  CML     &          -SINWAT*(GWATY(I,J,bi,bj)-VICEC(I,J,bi,bj)))             areaW = 0.5 * (AREA(I,J,bi,bj) + AREA(I-1,J,bi,bj))
65  CML           WINDY(I,J,bi,bj)=DWATN(I,J,bi,bj)       &          * SEAICEstressFactor
66  CML     &          *(SINWAT*(GWATX(I,J,bi,bj)-UICEC(I,J,bi,bj))             areaS = 0.5 * (AREA(I,J,bi,bj) + AREA(I,J-1,bi,bj))
67  CML     &          +COSWAT*(GWATY(I,J,bi,bj)-VICE(I,J,1,bi,bj)))       &          * SEAICEstressFactor
68  CML           WINDX(I,J,bi,bj)=WINDX(I,J,bi,bj)-( COR_ICE(I,J,bi,bj)             fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)
69  CML     &          *GWATY(I,J,bi,bj)-COR_ICE(I,J,bi,bj)*VICEC(I,J,bi,bj))       &          + areaW*taux(I,J,bi,bj)
70  CML           WINDY(I,J,bi,bj)=WINDY(I,J,bi,bj)-(-COR_ICE(I,J,bi,bj)       &          + stressDivergenceX(I,J,bi,bj) * SEAICEstressFactor
71  CML     &          *GWATX(I,J,bi,bj)+COR_ICE(I,J,bi,bj)*UICEC(I,J,bi,bj))             fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)
72  CML           WINDX(I,J,bi,bj)=WINDX(I,J,bi,bj)-(UICE(I,J,1,bi,bj)       &          + areaS*tauy(I,J,bi,bj)
73  CML     &          -UICE(I,J,3,bi,bj))*AMASS(I,J,bi,bj)/SEAICE_DT*TWO       &          + stressDivergenceY(I,J,bi,bj) * SEAICEstressFactor
74  CML           WINDY(I,J,bi,bj)=WINDY(I,J,bi,bj)-(VICE(I,J,1,bi,bj)            ENDDO
 CML     &          -VICE(I,J,3,bi,bj))*AMASS(I,J,bi,bj)/SEAICE_DT*TWO  
 CML          ENDDO  
 CML         ENDDO  
 CML        ENDDO  
 CML       ENDDO  
 CML       DO bj=myByLo(myThid),myByHi(myThid)  
 CML        DO bi=myBxLo(myThid),myBxHi(myThid)  
 CML         DO j=1,sNy  
 CML          DO i=1,sNx  
 CML           WINDX(I,J,bi,bj)=-WINDX(I,J,bi,bj)  
 CML           WINDY(I,J,bi,bj)=-WINDY(I,J,bi,bj)  
 CML          ENDDO  
 CML         ENDDO  
 CML        ENDDO  
 CML       ENDDO  
 CML      ENDIF  
 CML#endif /*  SEAICE_ALLOW_DYNAMICS */  
 CML#endif /* SEAICE_ORIGINAL_BAD_ICE_STRESS */  
   
 C--   Update overlap regions  
       CALL EXCH_UV_XY_RL(WINDX, WINDY, .TRUE., myThid)  
   
 #ifndef SEAICE_EXTERNAL_FLUXES  
 C--   Interpolate wind stress (N/m^2) from South-West B-grid  
 C     to South-West C-grid for forcing ocean model.  
       DO bj=myByLo(myThid),myByHi(myThid)  
        DO bi=myBxLo(myThid),myBxHi(myThid)  
         DO j=1,sNy  
          DO i=1,sNx  
             fu(I,J,bi,bj)=WINDX(I,J,bi,bj)  
             fv(I,J,bi,bj)=WINDY(I,J,bi,bj)  
75           ENDDO           ENDDO
76          ENDDO          ENDDO
77         ENDDO         ENDDO
       ENDDO  
       CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)  
 #endif /* ifndef SEAICE_EXTERNAL_FLUXES */  
78    
79  #ifdef SEAICE_TEST_ICE_STRESS_1        ELSE
80    C     else: useHB87StressCoupling=F
81    
82  C--   Compute ice-affected wind stress (interpolate to U/V-points)  C--   Compute ice-affected wind stress (interpolate to U/V-points)
83    C     by averaging wind stress and ice-ocean stress according to
84    C     ice cover
85        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
86         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
87          DO j=1,sNy          DO j=1,sNy
88           DO i=1,sNx           DO i=1,sNx
89            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) )*
90       &         COSWAT *       &         COSWAT *
91       &         ( UICE(I,J,1,bi,bj)-GWATX(I,J,bi,bj) )       &         ( uIce(I,J,bi,bj)-uVel(I,J,kSrf,bi,bj) )
92       &         - SINWAT* 0.5 _d 0 * (       &         - SIGN(SINWAT, _fCori(I,J,bi,bj)) * 0.5 _d 0 *
93       &          0.5 _d 0*(vIce(I  ,J  ,1,bi,bj)-GWATY(I  ,J  ,bi,bj)       &         ( DWATN(I  ,J,bi,bj) *
94       &                   +vIce(I-1,J  ,1,bi,bj)-GWATY(I-1,J  ,bi,bj))       &         0.5 _d 0*(vIce(I  ,J  ,bi,bj)-vVel(I  ,J  ,kSrf,bi,bj)
95       &         +0.5 _d 0*(vIce(I  ,J+1,1,bi,bj)-GWATY(I  ,J+1,bi,bj)       &                  +vIce(I  ,J+1,bi,bj)-vVel(I  ,J+1,kSrf,bi,bj))
96       &                   +vIce(I-1,J+1,1,bi,bj)-GWATY(I-1,J+1,bi,bj)) )       &         + DWATN(I-1,J,bi,bj) *
97         &         0.5 _d 0*(vIce(I-1,J  ,bi,bj)-vVel(I-1,J  ,kSrf,bi,bj)
98         &                  +vIce(I-1,J+1,bi,bj)-vVel(I-1,J+1,kSrf,bi,bj))
99       &         )       &         )
100            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) )*
101       &         SINWAT *       &         COSWAT *
102       &         ( UICE(I,J,1,bi,bj)-GWATX(I,J,bi,bj) )       &         ( vIce(I,J,bi,bj)-vVel(I,J,kSrf,bi,bj) )
103       &         + COSWAT * 0.5 _d 0 * (       &         + SIGN(SINWAT,  _fCori(I,J,bi,bj)) * 0.5 _d 0 *
104       &          0.5 _d 0*(uIce(I  ,J  ,1,bi,bj)-GWATX(I  ,J  ,bi,bj)       &         ( DWATN(I,J  ,bi,bj) *
105       &                   +uIce(I+1,J  ,1,bi,bj)-GWATX(I+1,J  ,bi,bj))       &         0.5 _d 0*(uIce(I  ,J  ,bi,bj)-uVel(I  ,J  ,kSrf,bi,bj)
106       &         +0.5 _d 0*(uIce(I  ,J-1,1,bi,bj)-GWATX(I  ,J-1,bi,bj)       &                  +uIce(I+1,J  ,bi,bj)-uVel(I+1,J  ,kSrf,bi,bj))
107       &                   +uIce(I+1,J-1,1,bi,bj)-GWATX(I+1,J-1,bi,bj)) )       &         + DWATN(I,J-1,bi,bj) *
108         &         0.5 _d 0*(uIce(I  ,J-1,bi,bj)-uVel(I  ,J-1,kSrf,bi,bj)
109         &                  +uIce(I+1,J-1,bi,bj)-uVel(I+1,J-1,kSrf,bi,bj))
110       &         )       &         )
111            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,bi,bj) + AREA(I-1,J,bi,bj))
112       &         AREA(I,J,1,bi,bj)*fuIce       &         * SEAICEstressFactor
113            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,bi,bj) + AREA(I,J-1,bi,bj))
114       &         AREA(I,J,1,bi,bj)*fvIce       &         * SEAICEstressFactor
115              fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)+areaW*fuIceLoc
116              fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)+areaS*fvIceLoc
117           ENDDO           ENDDO
118          ENDDO          ENDDO
119         ENDDO         ENDDO
120        ENDDO        ENDDO
121          ENDIF
122        CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)        CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)
123  #endif /* SEAICE_TEST_ICE_STRESS_1 */  
124  #endif /* not SEAICE_CGRID */  #endif /* SEAICE_CGRID */
125    
126        RETURN        RETURN
127        END        END
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.27
  ViewVC Help
Powered by ViewVC 1.1.22