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

Diff of /MITgcm/pkg/seaice/dynsolver.F

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

revision 1.15 by dimitri, Mon Dec 27 20:34:11 2004 UTC revision 1.16 by mlosch, Wed Feb 15 09:04:37 2006 UTC
# Line 17  C     === Global variables === Line 17  C     === Global variables ===
17  #include "SIZE.h"  #include "SIZE.h"
18  #include "EEPARAMS.h"  #include "EEPARAMS.h"
19  #include "PARAMS.h"  #include "PARAMS.h"
20    #include "GRID.h"
21  #include "FFIELDS.h"  #include "FFIELDS.h"
22  #include "SEAICE.h"  #include "SEAICE.h"
23  #include "SEAICE_GRID.h"  CML#include "SEAICE_GRID.h"
24  #include "SEAICE_PARAMS.h"  #include "SEAICE_PARAMS.h"
25  #include "SEAICE_FFIELDS.h"  #include "SEAICE_FFIELDS.h"
26    
# Line 62  C--   FIRST SET UP BASIC CONSTANTS Line 63  C--   FIRST SET UP BASIC CONSTANTS
63        RHOAIR=1.3 _d 0        RHOAIR=1.3 _d 0
64        ECCEN=TWO        ECCEN=TWO
65        ECM2=ONE/(ECCEN**2)        ECM2=ONE/(ECCEN**2)
66          RADIUS=rSphere
67        RADIUS=6370. _d 3        RADIUS=6370. _d 3
68        PSTAR=SEAICE_strength        PSTAR=SEAICE_strength
69    
# Line 82  C--   NOW SET UP MASS PER UNIT AREA AND Line 84  C--   NOW SET UP MASS PER UNIT AREA AND
84         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
85          DO j=1,sNy          DO j=1,sNy
86           DO i=1,sNx           DO i=1,sNx
87            AMASS(I,J,bi,bj)=RHOICE*QUART*(HEFF(i,j,1,bi,bj)            AMASS(I,J,bi,bj)=RHOICE*QUART*(
88       &           +HEFF(i-1,j,1,bi,bj)       &          HEFF(i,j  ,1,bi,bj) + HEFF(i-1,j  ,1,bi,bj)
89       &           +HEFF(i,j-1,1,bi,bj)       &         +HEFF(i,j-1,1,bi,bj) + HEFF(i-1,j-1,1,bi,bj) )
90       &           +HEFF(i-1,j-1,1,bi,bj))            COR_ICE(I,J,bi,bj)=AMASS(I,J,bi,bj) * _fCoriG(I,J,bi,bj)
           COR_ICE(I,J,bi,bj)=AMASS(I,J,bi,bj)  
      &         *TWO*OMEGA*SINEICE(I,J,bi,bj)  
91           ENDDO           ENDDO
92          ENDDO          ENDDO
93         ENDDO         ENDDO
# Line 185  C NOW DO PREDICTOR TIME STEP Line 185  C NOW DO PREDICTOR TIME STEP
185          DO j=1,sNy          DO j=1,sNy
186           DO i=1,sNx           DO i=1,sNx
187  C NOW EVALUATE STRAIN RATES  C NOW EVALUATE STRAIN RATES
188            E11(I,J,bi,bj)=HALF/(DXTICE(I,J,bi,bj)*CSTICE(I,J,bi,bj))            E11(I,J,bi,bj)=HALF* _recip_dxF(I,J,bi,bj) *
189       &                *(UICE(I+1,J+1,1,bi,bj)+UICE(I+1,J,1,bi,bj)       &                (UICE(I+1,J+1,1,bi,bj)+UICE(I+1,J,1,bi,bj)
190       &                -UICE(I,J+1,1,bi,bj)-UICE(I,J,1,bi,bj))       &                -UICE(I,  J+1,1,bi,bj)-UICE(I,  J,1,bi,bj))
191       &                -QUART*(VICE(I+1,J+1,1,bi,bj)+VICE(I,J+1,1,bi,bj)       &                -QUART*
192       &                +VICE(I,J,1,bi,bj)+VICE(I+1,J,1,bi,bj))       &                (VICE(I+1,J+1,1,bi,bj)+VICE(I,  J+1,1,bi,bj)
193       &                *TNGTICE(I,J,bi,bj)/RADIUS       &                +VICE(I,  J,  1,bi,bj)+VICE(I+1,J, 1,bi,bj))
194            E22(I,J,bi,bj)=HALF/DYTICE(I,J,bi,bj)  CML     &                * _tanPhiAtU(I,J,bi,bj)*recip_rSphere
195       &                *(VICE(I+1,J+1,1,bi,bj)+VICE(I,J+1,1,bi,bj)       &                * _tanPhiAtU(I,J,bi,bj)/RADIUS
196       &                -VICE(I+1,J,1,bi,bj)-VICE(I,J,1,bi,bj))            E22(I,J,bi,bj)=HALF* _recip_dyF(I,J,bi,bj) *
197            E12(I,J,bi,bj)=HALF*(HALF/DYTICE(I,J,bi,bj)       &                (VICE(I+1,J+1,1,bi,bj)+VICE(I,J+1,1,bi,bj)
198       &                *(UICE(I+1,J+1,1,bi,bj)+UICE(I,J+1,1,bi,bj)       &                -VICE(I+1,J,  1,bi,bj)-VICE(I,J,  1,bi,bj))
199       &                -UICE(I+1,J,1,bi,bj)-UICE(I,J,1,bi,bj))            E12(I,J,bi,bj)=HALF*(HALF * _recip_dyF(I,J,bi,bj) *
200       &                +HALF/(DXTICE(I,J,bi,bj)*CSTICE(I,J,bi,bj))       &                (UICE(I+1,J+1,1,bi,bj)+UICE(I,J+1,1,bi,bj)
201       &                *(VICE(I+1,J+1,1,bi,bj)+VICE(I+1,J,1,bi,bj)       &                -UICE(I+1,J,  1,bi,bj)-UICE(I,J,  1,bi,bj))
202       &                -VICE(I,J+1,1,bi,bj)-VICE(I,J,1,bi,bj))       &                +HALF * _recip_dxF(I,J,bi,bj) *
203       &                +QUART*(UICE(I+1,J+1,1,bi,bj)+UICE(I,J+1,1,bi,bj)       &                (VICE(I+1,J+1,1,bi,bj)+VICE(I+1,J,1,bi,bj)
204       &                +UICE(I,J,1,bi,bj)+UICE(I+1,J,1,bi,bj))       &                -VICE(I,  J+1,1,bi,bj)-VICE(I,  J,1,bi,bj))
205       &                *TNGTICE(I,J,bi,bj)/RADIUS)       &                +QUART*
206         &                (UICE(I+1,J+1,1,bi,bj)+UICE(I,  J+1,1,bi,bj)
207         &                +UICE(I,  J,  1,bi,bj)+UICE(I+1,J,  1,bi,bj))
208    CML     &                * _tanPhiAtU(I,J,bi,bj)*recip_rSphere)
209         &                * _tanPhiAtU(I,J,bi,bj)/RADIUS)
210  C  NOW EVALUATE VISCOSITIES  C  NOW EVALUATE VISCOSITIES
211            DELT1=(E11(I,J,bi,bj)**2+E22(I,J,bi,bj)**2)*(ONE+ECM2)            DELT1=(E11(I,J,bi,bj)**2+E22(I,J,bi,bj)**2)*(ONE+ECM2)
212       &        +4.0 _d 0*ECM2*E12(I,J,bi,bj)**2       &        +4.0 _d 0*ECM2*E12(I,J,bi,bj)**2
# Line 255  C NOW ADD IN CURRENT FORCE Line 259  C NOW ADD IN CURRENT FORCE
259       &                     +COSWAT*GWATY(I,J,bi,bj))       &                     +COSWAT*GWATY(I,J,bi,bj))
260  C NOW CALCULATE PRESSURE FORCE AND ADD TO EXTERNAL FORCE  C NOW CALCULATE PRESSURE FORCE AND ADD TO EXTERNAL FORCE
261            FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj)            FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj)
262       &          -(QUART/(DXUICE(I,J,bi,bj)*CSUICE(I,J,bi,bj)))       &          -QUART * _recip_dxV(I,J,bi,bj)
263       &          *(PRESS(I,J,bi,bj)+PRESS(I,J-1,bi,bj)       &          *(PRESS(I,  J,bi,bj) + PRESS(I,  J-1,bi,bj)
264       &          -PRESS(I-1,J,bi,bj)-PRESS(I-1,J-1,bi,bj))       &          - PRESS(I-1,J,bi,bj) - PRESS(I-1,J-1,bi,bj))
265            FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj)-QUART/DYUICE(I,J,bi,bj)            FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj)
266       &          *(PRESS(I,J,bi,bj)+PRESS(I-1,J,bi,bj)       &         -QUART * _recip_dyU(I,J,bi,bj)
267       &          -PRESS(I,J-1,bi,bj)-PRESS(I-1,J-1,bi,bj))       &          *(PRESS(I,J,  bi,bj) + PRESS(I-1,J,  bi,bj)
268         &          - PRESS(I,J-1,bi,bj) - PRESS(I-1,J-1,bi,bj))
269           ENDDO           ENDDO
270          ENDDO          ENDDO
271         ENDDO         ENDDO
# Line 292  C NOW DO MODIFIED EULER STEP Line 297  C NOW DO MODIFIED EULER STEP
297          DO j=1,sNy          DO j=1,sNy
298           DO i=1,sNx           DO i=1,sNx
299  C NOW EVALUATE STRAIN RATES  C NOW EVALUATE STRAIN RATES
300            E11(I,J,bi,bj)=HALF/(DXTICE(I,J,bi,bj)*CSTICE(I,J,bi,bj))            E11(I,J,bi,bj)=HALF * _recip_dxF(I,J,bi,bj) *
301       &                *(UICE(I+1,J+1,1,bi,bj)+UICE(I+1,J,1,bi,bj)       &                (UICE(I+1,J+1,1,bi,bj)+UICE(I+1,J,1,bi,bj)
302       &                -UICE(I,J+1,1,bi,bj)-UICE(I,J,1,bi,bj))       &                -UICE(I,  J+1,1,bi,bj)-UICE(I,  J,1,bi,bj))
303       &                -QUART*(VICE(I+1,J+1,1,bi,bj)+VICE(I,J+1,1,bi,bj)       &                -QUART*
304       &                +VICE(I,J,1,bi,bj)+VICE(I+1,J,1,bi,bj))       &                (VICE(I+1,J+1,1,bi,bj)+VICE(I,  J+1,1,bi,bj)
305       &                *TNGTICE(I,J,bi,bj)/RADIUS       &                +VICE(I,  J,  1,bi,bj)+VICE(I+1,J,  1,bi,bj))
306            E22(I,J,bi,bj)=HALF/DYTICE(I,J,bi,bj)  CML     &                * _tanPhiAtU(I,J,bi,bj)*recip_rSphere
307       &                *(VICE(I+1,J+1,1,bi,bj)+VICE(I,J+1,1,bi,bj)       &                * _tanPhiAtU(I,J,bi,bj)/RADIUS
308       &                -VICE(I+1,J,1,bi,bj)-VICE(I,J,1,bi,bj))            E22(I,J,bi,bj)=HALF * _recip_dyF(I,J,bi,bj) *
309            E12(I,J,bi,bj)=HALF*(HALF/DYTICE(I,J,bi,bj)       &                (VICE(I+1,J+1,1,bi,bj)+VICE(I,J+1,1,bi,bj)
310       &                *(UICE(I+1,J+1,1,bi,bj)+UICE(I,J+1,1,bi,bj)       &                -VICE(I+1,J,  1,bi,bj)-VICE(I,J,  1,bi,bj))
311       &                -UICE(I+1,J,1,bi,bj)-UICE(I,J,1,bi,bj))            E12(I,J,bi,bj)=HALF*(HALF * _recip_dyF(I,J,bi,bj) *
312       &                +HALF/(DXTICE(I,J,bi,bj)*CSTICE(I,J,bi,bj))       &                (UICE(I+1,J+1,1,bi,bj)+UICE(I,J+1,1,bi,bj)
313       &                *(VICE(I+1,J+1,1,bi,bj)+VICE(I+1,J,1,bi,bj)       &                -UICE(I+1,J,  1,bi,bj)-UICE(I,J,  1,bi,bj))
314       &                -VICE(I,J+1,1,bi,bj)-VICE(I,J,1,bi,bj))       &                +HALF * _recip_dxF(I,J,bi,bj) *
315       &                +QUART*(UICE(I+1,J+1,1,bi,bj)+UICE(I,J+1,1,bi,bj)       &                (VICE(I+1,J+1,1,bi,bj)+VICE(I+1,J,1,bi,bj)
316       &                +UICE(I,J,1,bi,bj)+UICE(I+1,J,1,bi,bj))       &                -VICE(I,  J+1,1,bi,bj)-VICE(I,  J,1,bi,bj))
317       &                *TNGTICE(I,J,bi,bj)/RADIUS)       &                +QUART*
318         &                (UICE(I+1,J+1,1,bi,bj)+UICE(I,  J+1,1,bi,bj)
319         &                +UICE(I,  J,  1,bi,bj)+UICE(I+1,J,  1,bi,bj))
320    CML     &                * _tanPhiAtU(I,J,bi,bj)*recip_rSphere)
321         &                * _tanPhiAtU(I,J,bi,bj)/RADIUS)
322  C  NOW EVALUATE VISCOSITIES  C  NOW EVALUATE VISCOSITIES
323            DELT1=(E11(I,J,bi,bj)**2+E22(I,J,bi,bj)**2)*(ONE+ECM2)            DELT1=(E11(I,J,bi,bj)**2+E22(I,J,bi,bj)**2)*(ONE+ECM2)
324       &        +4. _d 0*ECM2*E12(I,J,bi,bj)**2       &        +4. _d 0*ECM2*E12(I,J,bi,bj)**2
# Line 362  C NOW ADD IN CURRENT FORCE Line 371  C NOW ADD IN CURRENT FORCE
371       &                     +COSWAT*GWATY(I,J,bi,bj))       &                     +COSWAT*GWATY(I,J,bi,bj))
372  C NOW CALCULATE PRESSURE FORCE AND ADD TO EXTERNAL FORCE  C NOW CALCULATE PRESSURE FORCE AND ADD TO EXTERNAL FORCE
373            FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj)            FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj)
374       &          -(QUART/(DXUICE(I,J,bi,bj)*CSUICE(I,J,bi,bj)))       &          -QUART * _recip_dxV(I,J,bi,bj)
375       &          *(PRESS(I,J,bi,bj)+PRESS(I,J-1,bi,bj)       &          *(PRESS(I,  J,bi,bj) + PRESS(I,  J-1,bi,bj)
376       &          -PRESS(I-1,J,bi,bj)-PRESS(I-1,J-1,bi,bj))       &          - PRESS(I-1,J,bi,bj) - PRESS(I-1,J-1,bi,bj))
377            FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj)-QUART/DYUICE(I,J,bi,bj)            FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj)
378       &          *(PRESS(I,J,bi,bj)+PRESS(I-1,J,bi,bj)       &         -QUART * _recip_dyU(I,J,bi,bj)
379       &          -PRESS(I,J-1,bi,bj)-PRESS(I-1,J-1,bi,bj))       &          *(PRESS(I,J,  bi,bj) + PRESS(I-1,J,  bi,bj)
380         &          - PRESS(I,J-1,bi,bj) - PRESS(I-1,J-1,bi,bj))
381           ENDDO           ENDDO
382          ENDDO          ENDDO
383         ENDDO         ENDDO
Legend:
Removed from v.1.15  
changed lines
  Added in v.1.16
  ViewVC Help
Powered by ViewVC 1.1.22