/[MITgcm]/MITgcm/model/src/ini_cg3d.F
ViewVC logotype

Diff of /MITgcm/model/src/ini_cg3d.F

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

revision 1.18 by jmc, Tue Nov 8 06:06:07 2005 UTC revision 1.19 by jmc, Tue Dec 20 20:31:28 2005 UTC
# Line 25  C     === Global variables === Line 25  C     === Global variables ===
25  #include "EEPARAMS.h"  #include "EEPARAMS.h"
26  #include "PARAMS.h"  #include "PARAMS.h"
27  #include "GRID.h"  #include "GRID.h"
28    #include "SURFACE.h"
29  #include "CG3D.h"  #include "CG3D.h"
30  #include "SOLVE_FOR_PRESSURE3D.h"  #include "SOLVE_FOR_PRESSURE3D.h"
31  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
# Line 46  C     faceArea - Temporary used to hold Line 47  C     faceArea - Temporary used to hold
47  C     myNorm - Work variable used in clculating normalisation factor  C     myNorm - Work variable used in clculating normalisation factor
48        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
49        INTEGER bi, bj        INTEGER bi, bj
50        INTEGER I,  J, K        INTEGER I,  J, K, ks
51        _RL     faceArea        _RL     faceArea
52        _RS     myNorm        _RS     myNorm
53        _RL     theRecip_Dr        _RL     theRecip_Dr
# Line 67  C-       From common bloc CG3D_R: Line 68  C-       From common bloc CG3D_R:
68             aW3d(I,J,K,bi,bj) = 0.             aW3d(I,J,K,bi,bj) = 0.
69             aS3d(I,J,K,bi,bj) = 0.             aS3d(I,J,K,bi,bj) = 0.
70             aV3d(I,J,K,bi,bj) = 0.             aV3d(I,J,K,bi,bj) = 0.
71               aC3d(I,J,K,bi,bj) = 0.
72             zMC (I,J,K,bi,bj) = 0.             zMC (I,J,K,bi,bj) = 0.
73             zML (I,J,K,bi,bj) = 0.             zML (I,J,K,bi,bj) = 0.
74             zMU (I,J,K,bi,bj) = 0.             zMU (I,J,K,bi,bj) = 0.
# Line 83  C-       From common bloc CG3D_WK_R: Line 85  C-       From common bloc CG3D_WK_R:
85  C-       From common bloc SFP3D_COMMON_R:  C-       From common bloc SFP3D_COMMON_R:
86           DO J=1-OLy,sNy+OLy           DO J=1-OLy,sNy+OLy
87            DO I=1-OLx,sNx+OLx            DO I=1-OLx,sNx+OLx
 c          cg3d_x(I,J,K,bi,bj) = 0.  
88             cg3d_b(I,J,K,bi,bj) = 0.             cg3d_b(I,J,K,bi,bj) = 0.
89            ENDDO            ENDDO
90           ENDDO           ENDDO
# Line 100  C     aV3d: Ar/dR Line 101  C     aV3d: Ar/dR
101         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
102          DO K=1,Nr          DO K=1,Nr
103           DO J=1,sNy           DO J=1,sNy
104            DO I=1,sNx            DO I=1,sNx+1
105             faceArea = _dyG(I,J,bi,bj)*drF(K)             faceArea = _dyG(I,J,bi,bj)*drF(K)
106       &                *_hFacW(I,J,K,bi,bj)       &                *_hFacW(I,J,K,bi,bj)
107             aW3d(I,J,K,bi,bj) = faceArea*recip_dxC(I,J,bi,bj)             aW3d(I,J,K,bi,bj) = faceArea*recip_dxC(I,J,bi,bj)
108               myNorm = MAX(ABS(aW3d(I,J,K,bi,bj)),myNorm)
109              ENDDO
110             ENDDO
111             DO J=1,sNy+1
112              DO I=1,sNx
113             faceArea = _dxG(I,J,bi,bj)*drF(K)             faceArea = _dxG(I,J,bi,bj)*drF(K)
114       &                *_hFacS(I,J,K,bi,bj)       &                *_hFacS(I,J,K,bi,bj)
115             aS3d(I,J,K,bi,bj) = faceArea*recip_dyC(I,J,bi,bj)             aS3d(I,J,K,bi,bj) = faceArea*recip_dyC(I,J,bi,bj)
            myNorm = MAX(ABS(aW3d(I,J,K,bi,bj)),myNorm)  
116             myNorm = MAX(ABS(aS3d(I,J,K,bi,bj)),myNorm)             myNorm = MAX(ABS(aS3d(I,J,K,bi,bj)),myNorm)
117            ENDDO            ENDDO
118           ENDDO           ENDDO
# Line 195  CcnhDebugEnds Line 200  CcnhDebugEnds
200        _END_MASTER( myThid )        _END_MASTER( myThid )
201        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
202         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
203    C-    Set solver main diagonal term
204            DO K=1,Nr
205             DO J=1,sNy
206              DO I=1,sNx
207               aW = aW3d(I  ,J,K,bi,bj)
208               aE = aW3d(I+1,J,K,bi,bj)
209               aN = aS3d(I,J+1,K,bi,bj)
210               aS = aS3d(I,J  ,K,bi,bj)
211               aU = aV3d(I,J,K,bi,bj)
212               IF ( K .NE. Nr  ) THEN
213                aL = aV3d(I,J,K+1,bi,bj)
214               ELSE
215                aL = 0.
216               ENDIF
217               aC3d(I,J,K,bi,bj) = -aW-aE-aN-aS-aU-aL
218              ENDDO
219             ENDDO
220            ENDDO
221    C-    Add free-surface source term
222            DO J=1,sNy
223             DO I=1,sNx
224               ks = ksurfC(I,J,bi,bj)
225               IF ( ks.LE.Nr ) THEN
226                 aC3d(I,J,ks,bi,bj) = aC3d(I,J,ks,bi,bj)
227         &         -freeSurfFac*recip_Bo(I,J,bi,bj)
228         &          *rA(I,J,bi,bj)/deltaTMom/deltaTfreesurf
229               ENDIF
230             ENDDO
231            ENDDO
232    C-    Matrix solver normalisation
233          DO K=1,Nr          DO K=1,Nr
234           DO J=1,sNy           DO J=1,sNy
235            DO I=1,sNx            DO I=1,sNx
236             aW3d(I,J,K,bi,bj) = aW3d(I,J,K,bi,bj)*myNorm             aW3d(I,J,K,bi,bj) = aW3d(I,J,K,bi,bj)*myNorm
237             aS3d(I,J,K,bi,bj) = aS3d(I,J,K,bi,bj)*myNorm             aS3d(I,J,K,bi,bj) = aS3d(I,J,K,bi,bj)*myNorm
238             aV3d(I,J,K,bi,bj) = aV3d(I,J,K,bi,bj)*myNorm             aV3d(I,J,K,bi,bj) = aV3d(I,J,K,bi,bj)*myNorm
239               aC3d(I,J,K,bi,bj) = aC3d(I,J,K,bi,bj)*myNorm
240            ENDDO            ENDDO
241           ENDDO           ENDDO
242          ENDDO          ENDDO
243         ENDDO         ENDDO
244        ENDDO        ENDDO
245          
246  C--   Update overlap regions  C--   Update overlap regions
247  c     _EXCH_XYZ_R4(aW3d, myThid)  c     _EXCH_XYZ_R4(aW3d, myThid)
248  c     _EXCH_XYZ_R4(aS3d, myThid)  c     _EXCH_XYZ_R4(aS3d, myThid)
249        CALL EXCH_UV_XYZ_RS(aW3d,aS3d,.FALSE.,myThid)        CALL EXCH_UV_XYZ_RS(aW3d,aS3d,.FALSE.,myThid)
250        _EXCH_XYZ_R4(aV3d, myThid)        _EXCH_XYZ_R4(aV3d, myThid)
251          _EXCH_XYZ_R4(aC3d, myThid)
252  CcnhDebugStarts  CcnhDebugStarts
253  C     CALL PLOT_FIELD_XYZRS( aW3d, 'AW3D INI_CG3D.1' , Nr, 1, myThid )  C     CALL PLOT_FIELD_XYZRS( aW3d, 'AW3D INI_CG3D.1' , Nr, 1, myThid )
254  C     CALL PLOT_FIELD_XYZRS( aS3d, 'AS3D INI_CG3D.1' , Nr, 1, myThid )  C     CALL PLOT_FIELD_XYZRS( aS3d, 'AS3D INI_CG3D.1' , Nr, 1, myThid )
# Line 224  C     check for consistency with S/R CG3 Line 261  C     check for consistency with S/R CG3
261  C     assuming zML is lower and zMU is upper!  C     assuming zML is lower and zMU is upper!
262        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
263         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
264          DO K=1,Nr            DO K=1,Nr
265           DO J=1,sNy           DO J=1,sNy
266            DO I=1,sNx            DO I=1,sNx
267             aW = aW3d(I  ,J,K,bi,bj)             IF ( aC3d(I,J,K,bi,bj) .NE. 0. ) THEN
268             aE = aW3d(I+1,J,K,bi,bj)              zMC(i,j,k,bi,bj) = aC3d(I,J,K,bi,bj)
269             aN = aS3d(I,J+1,K,bi,bj)              zML(i,j,k,bi,bj) = aV3d(I,J,K,bi,bj)
270             aS = aS3d(I,J  ,K,bi,bj)              zMU(i,j,k,bi,bj) = 0.
271             IF ( K .NE. 1 ) THEN              IF ( K.NE.Nr )
272              aU = aV3d(I,J,K,bi,bj)       &      zMU(i,j,k,bi,bj) = aV3d(I,J,K+1,bi,bj)
            ELSE  
             aU = 0  
            ENDIF  
            IF ( K .NE. Nr ) THEN  
             aL = aV3d(I,J,K+1,bi,bj)  
            ELSE  
             aL = 0  
            ENDIF  
            aC = -aW-aE-aN-aS-aU-aL  
            IF ( K .EQ. 1 .AND. aC .NE. 0. ) THEN  
             aC = aC  
      &         -freeSurfFac*myNorm*(horiVertRatio/gravity)*  
      &          rA(I,J,bi,bj)/deltaTMom/deltaTMom  
            ENDIF  
            IF ( aC .NE. 0. ) THEN  
             zMC(i,j,k,bi,bj) = aC  
             zMU(i,j,k,bi,bj) = aL  
             zML(i,j,k,bi,bj) = aU  
273  CcnhDebugStarts  CcnhDebugStarts
274  C           zMC(i,j,k,bi,bj) = 1.  C           zMC(i,j,k,bi,bj) = 1.
275  C           zMU(i,j,k,bi,bj) = 0.  C           zMU(i,j,k,bi,bj) = 0.
# Line 272  CcnhDebugEnds Line 291  CcnhDebugEnds
291       &     zMU(i,j,1,bi,bj)*zMC(i,j,1,bi,bj)       &     zMU(i,j,1,bi,bj)*zMC(i,j,1,bi,bj)
292           ENDDO           ENDDO
293          ENDDO          ENDDO
294          DO K=2,Nr            DO K=2,Nr
295           DO J=1,sNy           DO J=1,sNy
296            DO I=1,sNx            DO I=1,sNx
297             zMC(i,j,k,bi,bj) = 1. _d 0 /             zMC(i,j,k,bi,bj) = 1. _d 0 /
# Line 281  CcnhDebugEnds Line 300  CcnhDebugEnds
300            ENDDO            ENDDO
301           ENDDO           ENDDO
302          ENDDO          ENDDO
303          DO K=1,Nr            DO K=1,Nr
304           DO J=1,sNy           DO J=1,sNy
305            DO I=1,sNx            DO I=1,sNx
306             aW = aW3d(I  ,J,K,bi,bj)             aW = aW3d(I  ,J,K,bi,bj)
# Line 291  CcnhDebugEnds Line 310  CcnhDebugEnds
310             IF ( K .NE. 1 ) THEN             IF ( K .NE. 1 ) THEN
311              aU = aV3d(I,J,K-1,bi,bj)              aU = aV3d(I,J,K-1,bi,bj)
312             ELSE             ELSE
313              aU = 0              aU = 0.
314             ENDIF             ENDIF
315             IF ( K .NE. Nr ) THEN             IF ( K .NE. Nr ) THEN
316              aL = aV3d(I,J,K+1,bi,bj)              aL = aV3d(I,J,K+1,bi,bj)
317             ELSE             ELSE
318              aL = 0              aL = 0.
319             ENDIF             ENDIF
320             aC = -aW-aE-aN-aS-aU-aL             aC = -aW-aE-aN-aS-aU-aL
321             IF ( aC .EQ. 0. ) THEN             IF ( aC .EQ. 0. ) THEN
Legend:
Removed from v.1.18  
changed lines
  Added in v.1.19
  ViewVC Help
Powered by ViewVC 1.1.22