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

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

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

revision 1.4 by mlosch, Mon Jun 5 22:50:37 2006 UTC revision 1.5 by jmc, Wed Nov 1 01:56:23 2006 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "SEAICE_OPTIONS.h"  #include "SEAICE_OPTIONS.h"
5    
6  CStartOfInterface  CBOP
7        SUBROUTINE SEAICE_DIFFUSION(  C     !ROUTINE: SEAICE_DIFFUSION
8       U     HEFF,  C     !INTERFACE:
9       I     HEFFM, DELTT, myTime, myIter, myThid )        SUBROUTINE SEAICE_DIFFUSION(
10  C     /==========================================================\       I     tracerIdentity,
11  C     | SUBROUTINE advect                                        |       I     iceFld, iceMask, xA, yA,
12  C     | o Calculate ice advection                                |       U     gFld,
13  C     |==========================================================|       I     bi, bj, myTime, myIter, myThid )
14  C     \==========================================================/  
15    C     !DESCRIPTION: \bv
16    C     *==========================================================*
17    C     | SUBROUTINE SEAICE_DIFFUSION
18    C     | o Add tendency from horizontal diffusion
19    C     *==========================================================*
20    C     *==========================================================*
21    
22    C     !USES:
23        IMPLICIT NONE        IMPLICIT NONE
24    
25  C     === Global variables ===  C     === Global variables ===
# Line 26  CML#include "SEAICE_GRID.h" Line 34  CML#include "SEAICE_GRID.h"
34  # include "tamc.h"  # include "tamc.h"
35  #endif  #endif
36    
37    C     !INPUT PARAMETERS:
38  C     === Routine arguments ===  C     === Routine arguments ===
39  C     myThid - Thread no. that called this routine.  C     afx        :: horizontal advective flux, x direction
40        _RL HEFF       (1-OLx:sNx+OLx,1-OLy:sNy+OLy,3,nSx,nSy)  C     afy        :: horizontal advective flux, y direction
41        _RL HEFFM      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,  nSx,nSy)  C     myThid     :: my Thread Id number
42        _RL     DELTT        _RL iceFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
43          _RL gFld   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
44          _RL iceMask(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
45          _RS xA     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
46          _RS yA     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
47          INTEGER advScheme
48          INTEGER tracerIdentity
49          INTEGER bi,bj
50        _RL     myTime        _RL     myTime
51        INTEGER myIter        INTEGER myIter
52        INTEGER myThid        INTEGER myThid
53  CEndOfInterface  CEOP
54    
55    C     !LOCAL VARIABLES:
56  C     === Local variables ===  C     === Local variables ===
57  C     i,j,k,bi,bj - Loop counters  C     i,j :: Loop counters
58          INTEGER i, j, k
59        INTEGER i, j, bi, bj        _RL fZon   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
60          _RL fMer   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
61        _RL DIFFA(1-OLx:sNx+OLx, 1-OLy:sNy+OLy,nSx,nSy)  
62          IF ( diff1 .GT. 0. _d 0 ) THEN
63  #ifdef ALLOW_AUTODIFF_TAMC  C--   Tendency due to horizontal diffusion
64  cphCADJ STORE heff  = comlev1, key = ikey_dynamics          k = 1
65  #endif /* ALLOW_AUTODIFF_TAMC */          DO j=1-Oly,sNy+Oly
66             DO i=1-Olx,sNx+Olx
67              fZon  (i,j) = 0. _d 0
68  C--   This would be the natural way to do diffusion (explicitly)            fMer  (i,j) = 0. _d 0
 C     For now we stick to the modified Eulerian time step  
 CML      DO bj=myByLo(myThid),myByHi(myThid)  
 CML       DO bi=myBxLo(myThid),myBxHi(myThid)  
 CML        CALL GAD_DIFF_X(bi,bj,k,xA,diff1,localT,df,myThid)  
 CML        DO j=1-Oly,sNy+Oly  
 CML         DO i=1-Olx,sNx+Olx  
 CML          fZon(i,j) = fZon(i,j) + df(i,j)  
 CML         ENDDO  
 CML        ENDDO  
 CML        CALL GAD_DIFF_Y(bi,bj,k,yA,diff1,localT,df,myThid)  
 CML        DO j=1-Oly,sNy+Oly  
 CML         DO i=1-Olx,sNx+Olx  
 CML          fMer(i,j) = fMer(i,j) + df(i,j)  
 CML         ENDDO  
 CML        ENDDO  
 CMLC--   Divergence of fluxes: update scalar field  
 CML        DO j=1-Oly,sNy+Oly-1  
 CML         DO i=1-Olx,sNx+Olx-1  
 CML          HEFF(i,j,1,bi,bj)=HEFF(i,j,1,bi,bj) + DELTT *  
 CML     &         maskC(i,j,kSurface,bi,bj)*recip_rA(i,j,bi,bj)  
 CML     &       *( (fZon(i+1,j)-fZon(i,j))  
 CML     &         +(fMer(i,j+1)-fMer(i,j))  
 CML     &                       )  
 CML     &         )  
 CML         ENDDO  
 CML        ENDDO  
 CML       ENDDO  
 CML      ENDDO  
   
 C NOW DO DIFFUSION ON H(I,J,3)  
 C NOW CALCULATE DIFFUSION COEF ROUGHLY  
       DO bj=myByLo(myThid),myByHi(myThid)  
        DO bi=myBxLo(myThid),myBxHi(myThid)  
         DO j=1-OLy,sNy+OLy  
          DO i=1-OLx,sNx+OLx  
           DIFFA(I,J,bi,bj)=  
      &         DIFF1*MIN( _dxF(I,J,bi,bj), _dyF(I,J,bi,bj))  
          ENDDO  
         ENDDO  
        ENDDO  
       ENDDO  
       CALL DIFFUS(HEFF,DIFFA,HEFFM,DELTT, myThid)  
   
       DO bj=myByLo(myThid),myByHi(myThid)  
        DO bi=myBxLo(myThid),myBxHi(myThid)  
         DO j=1-OLy,sNy+OLy  
          DO i=1-OLx,sNx+OLx  
           HEFF(I,J,1,bi,bj)=(HEFF(I,J,1,bi,bj)+HEFF(I,J,3,bi,bj))  
      &                      *HEFFM(I,J,bi,bj)  
          ENDDO  
         ENDDO  
        ENDDO  
       ENDDO  
   
 C NOW CALCULATE DIFFUSION COEF ROUGHLY  
       DO bj=myByLo(myThid),myByHi(myThid)  
        DO bi=myBxLo(myThid),myBxHi(myThid)  
         DO j=1-OLy,sNy+OLy  
          DO i=1-OLx,sNx+OLx  
           DIFFA(I,J,bi,bj)=  
      &         -(MIN( _dxF(I,J,bi,bj),  _dyF(I,J,bi,bj)))**2/DELTT  
69           ENDDO           ENDDO
70          ENDDO          ENDDO
71         ENDDO  C--   X-direction
72        ENDDO          CALL GAD_DIFF_X(bi,bj,k,xA,diff1,iceFld,fZon,myThid)
73        CALL DIFFUS(HEFF,DIFFA,HEFFM,DELTT, myThid)  C--   Y-direction
74            CALL GAD_DIFF_Y(bi,bj,k,yA,diff1,iceFld,fMer,myThid)
75        DO bj=myByLo(myThid),myByHi(myThid)  C--   Divergence of fluxes: update scalar field
76         DO bi=myBxLo(myThid),myBxHi(myThid)  C--   Ugly:
77          DO j=1-OLy,sNy+OLy  C--   Apply factor min(DX,DY) to effectively end up with approximately
78           DO i=1-OLx,sNx+OLx  C--   the same diffusion coefficient as in subroutine ADVECT.
79            HEFF(I,J,1,bi,bj)=(HEFF(I,J,1,bi,bj)+HEFF(I,J,3,bi,bj))  C--   One day, I would like to rewrite the second order central difference
80       &                      *HEFFM(I,J,bi,bj)  C--   part, too, so that the value of DIFF1 has the same meaning as,
81    C--   say, diffKhT
82            DO j=1-Oly,sNy+Oly-1
83             DO i=1-Olx,sNx+Olx-1
84              gFld(i,j)= gFld(i,j)
85         &        - iceMask(i,j,bi,bj)*recip_rA(i,j,bi,bj)
86         &        *( (fZon(i+1,j)-fZon(i,j))
87         &         + (fMer(i,j+1)-fMer(i,j)) )
88         &        *MIN( _dxF(I,J,bi,bj), _dyF(I,J,bi,bj))
89           ENDDO           ENDDO
90          ENDDO          ENDDO
91         ENDDO  C     endif do horizontal diffusion
92        ENDDO        ENDIF
93    
94        RETURN        RETURN
95        END        END
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.22