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

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

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

revision 1.18 by cnh, Sun Sep 6 14:45:11 1998 UTC revision 1.31 by jmc, Tue Mar 6 16:51:02 2001 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    
# Line 11  C     |================================= Line 12  C     |=================================
12  C     | These arrays are purely a function of the basin geom.    |  C     | These arrays are purely a function of the basin geom.    |
13  C     | We set then here once and them use then repeatedly.      |  C     | We set then here once and them use then repeatedly.      |
14  C     \==========================================================/  C     \==========================================================/
15          IMPLICIT NONE
16    
17  C     === Global variables ===  C     === Global variables ===
18  #include "SIZE.h"  #include "SIZE.h"
# Line 18  C     === Global variables === Line 20  C     === Global variables ===
20  #include "PARAMS.h"  #include "PARAMS.h"
21  #include "GRID.h"  #include "GRID.h"
22  #include "DYNVARS.h"  #include "DYNVARS.h"
23  #include "CG2D.h"  #include "SURFACE.h"
24    #include "CG2D_INTERNAL.h"
25    #ifdef ALLOW_OBCS
26    #include "OBCS.h"
27    #endif
28    
29  C     === Routine arguments ===  C     === Routine arguments ===
30  C     myThid - Thread no. that called this routine.  C     myThid - Thread no. that called this routine.
# Line 40  C     myNorm - Work variable used in clc Line 46  C     myNorm - Work variable used in clc
46        _RS     myNorm        _RS     myNorm
47        _RL     aC, aCw, aCs        _RL     aC, aCw, aCs
48    
49    C-- Initialise -Boyancy at surface level : Bo_surf
50    C   Bo_surf is defined as d/dr(Phi_surf) and set to g/rtoz (linear free surface)
51    C     with rtoz = conversion factor from r-unit to z-unit  (=horiVertRatio)
52    C   an acurate formulation will include P_surf and T,S_surf effects on rho_surf:
53    C    z-ocean (rtoz=1) : Bo_surf = - Boyancy = gravity * rho_surf/rho_0
54    C    p-atmos (rtoz=rho_c*g) : Bo_surf = (1/rho)_surf
55    C--
56          IF ( buoyancyRelation .EQ. 'ATMOSPHERIC' ) THEN
57            DO bj=myByLo(myThid),myByHi(myThid)
58             DO bi=myBxLo(myThid),myBxHi(myThid)
59              DO J=1-Oly,sNy+Oly
60               DO I=1-Olx,sNx+Olx
61                 Bo_surf(I,J,bi,bj) = recip_rhoConst
62                 recip_Bo(I,J,bi,bj) = rhoConst
63               ENDDO
64              ENDDO
65             ENDDO
66            ENDDO
67          ELSE
68            DO bj=myByLo(myThid),myByHi(myThid)
69             DO bi=myBxLo(myThid),myBxHi(myThid)
70              DO J=1-Oly,sNy+Oly
71               DO I=1-Olx,sNx+Olx
72                 Bo_surf(I,J,bi,bj) = gravity
73                 recip_Bo(I,J,bi,bj) = recip_Gravity
74               ENDDO
75              ENDDO
76             ENDDO
77            ENDDO
78          ENDIF
79    
80    C--   Update overlap regions
81          _EXCH_XY_R8(Bo_surf, myThid)
82          _EXCH_XY_R8(recip_Bo, myThid)
83    
84  C--   Initialise laplace operator  C--   Initialise laplace operator
85  C     aW2d: integral in Z Ax/dX  C     aW2d: integral in Z Ax/dX
86  C     aS2d: integral in Z Ay/dY  C     aS2d: integral in Z Ay/dY
# Line 55  C     aS2d: integral in Z Ay/dY Line 96  C     aS2d: integral in Z Ay/dY
96          DO K=1,Nr          DO K=1,Nr
97           DO J=1,sNy           DO J=1,sNy
98            DO I=1,sNx            DO I=1,sNx
99             faceArea = _dyG(I,J,bi,bj)*drF(K)*_hFacW(I,J,K,bi,bj)             faceArea = _dyG(I,J,bi,bj)*drF(K)
100         &               *_hFacW(I,J,K,bi,bj)
101             aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) +             aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) +
102         &      implicSurfPress*implicDiv2DFlow*
103       &      gBaro*faceArea*recip_dxC(I,J,bi,bj)       &      gBaro*faceArea*recip_dxC(I,J,bi,bj)
104             faceArea = _dxG(I,J,bi,bj)*drF(K)*_hFacS(I,J,K,bi,bj)             faceArea = _dxG(I,J,bi,bj)*drF(K)
105         &               *_hFacS(I,J,K,bi,bj)
106             aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) +             aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) +
107         &      implicSurfPress*implicDiv2DFlow*
108       &      gBaro*faceArea*recip_dyC(I,J,bi,bj)       &      gBaro*faceArea*recip_dyC(I,J,bi,bj)
109            ENDDO            ENDDO
110           ENDDO           ENDDO
111          ENDDO          ENDDO
112    #ifdef ALLOW_OBCS
113            IF (useOBCS) THEN
114             DO I=1,sNx
115              IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj),bi,bj)=0.
116              IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj)+1,bi,bj)=0.
117              IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj)+1,bi,bj)=0.
118              IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj),bi,bj)=0.
119             ENDDO
120             DO J=1,sNy
121              IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj),J,bi,bj)=0.
122              IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj)+1,J,bi,bj)=0.
123              IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj)+1,J,bi,bj)=0.
124              IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj),J,bi,bj)=0.
125             ENDDO
126            ENDIF
127    #endif
128          DO J=1,sNy          DO J=1,sNy
129           DO I=1,sNx           DO I=1,sNx
130            myNorm = MAX(ABS(aW2d(I,J,bi,bj)),myNorm)            myNorm = MAX(ABS(aW2d(I,J,bi,bj)),myNorm)
# Line 72  C     aS2d: integral in Z Ay/dY Line 133  C     aS2d: integral in Z Ay/dY
133          ENDDO          ENDDO
134         ENDDO         ENDDO
135        ENDDO        ENDDO
136        cg2dNbuf(1,myThid) = myNorm        _GLOBAL_MAX_R4( myNorm, myThid )
137        _GLOBAL_MAX_R4( cg2dNbuf, myNorm, myThid )        IF ( myNorm .NE. 0. _d 0 ) THEN
138        IF ( cg2dNbuf(1,1) .NE. 0. _d 0 ) THEN         myNorm = 1. _d 0/myNorm
        myNorm = 1. _d 0/cg2dNbuf(1,1)  
139        ELSE        ELSE
140         myNorm = 1. _d 0         myNorm = 1. _d 0
141        ENDIF        ENDIF
142         cg2dNorm = myNorm         cg2dNorm = myNorm
143        _BEGIN_MASTER( myThid )        _BEGIN_MASTER( myThid )
144  CcnhDebugStarts  CcnhDebugStarts
145         WRITE(msgBuf,'(A,F)') '// CG2D normalisation factor = ', cg2dNorm         WRITE(msgBuf,'(A,E40.25)') '// CG2D normalisation factor = ',
146         &               cg2dNorm
147         CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)         CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
148         WRITE(msgBuf,*) '                               '         WRITE(msgBuf,*) '                               '
149         CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)         CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
# Line 107  CcnhDebugEnds Line 168  CcnhDebugEnds
168        _EXCH_XY_R4(aW2d, myThid)        _EXCH_XY_R4(aW2d, myThid)
169        _EXCH_XY_R4(aS2d, myThid)        _EXCH_XY_R4(aS2d, myThid)
170  CcnhDebugStarts  CcnhDebugStarts
171        CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )
172        CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.2' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.2' , 1, myThid )
173  CcnhDebugEnds  CcnhDebugEnds
174    
175  C--   Initialise preconditioner  C--   Initialise preconditioner
# Line 132  C           defaults to 0.51 but can be Line 193  C           defaults to 0.51 but can be
193            pC(I,J,bi,bj) = 1. _d 0            pC(I,J,bi,bj) = 1. _d 0
194            aC = -(            aC = -(
195       &     aW2d(I,J,bi,bj) + aW2d(I+1,J  ,bi,bj)       &     aW2d(I,J,bi,bj) + aW2d(I+1,J  ,bi,bj)
196       &    +aS2d(I,J,bi,bj) + aS2D(I  ,J+1,bi,bj)       &    +aS2d(I,J,bi,bj) + aS2d(I  ,J+1,bi,bj)
197       &    +freeSurfFac*myNorm*  Gravity*rhoConst*       &    +freeSurfFac*myNorm* horiVertRatio*
198       &     rA(I,J,bi,bj)/deltaTMom/deltaTMom       &     rA(I,J,bi,bj)/deltaTMom/deltaTMom
199       &    )       &    )
200            aCs = -(            aCs = -(
201       &     aW2d(I,J-1,bi,bj) + aW2d(I+1,J-1,bi,bj)       &     aW2d(I,J-1,bi,bj) + aW2d(I+1,J-1,bi,bj)
202       &    +aS2d(I,J-1,bi,bj) + aS2d(I  ,J  ,bi,bj)       &    +aS2d(I,J-1,bi,bj) + aS2d(I  ,J  ,bi,bj)
203       &    +freeSurfFac*myNorm*  Gravity*rhoConst*       &    +freeSurfFac*myNorm* horiVertRatio*
204       &     rA(I,J-1,bi,bj)/deltaTMom/deltaTMom       &     rA(I,J-1,bi,bj)/deltaTMom/deltaTMom
205       &    )       &    )
206            aCw = -(            aCw = -(
207       &     aW2d(I-1,J,bi,bj) + aW2d(I  ,J  ,bi,bj)       &     aW2d(I-1,J,bi,bj) + aW2d(I  ,J  ,bi,bj)
208       &    +aS2d(I-1,J,bi,bj) + aS2d(I-1,J+1,bi,bj)       &    +aS2d(I-1,J,bi,bj) + aS2d(I-1,J+1,bi,bj)
209       &    +freeSurfFac*myNorm*  Gravity*rhoConst*       &    +freeSurfFac*myNorm* horiVertRatio*
210       &     rA(I-1,J,bi,bj)/deltaTMom/deltaTMom       &     rA(I-1,J,bi,bj)/deltaTMom/deltaTMom
211       &    )       &    )
212            IF ( aC .EQ. 0. ) THEN            IF ( aC .EQ. 0. ) THEN
213              pC(I,J,bi,bj) = 0. _d 0              pC(I,J,bi,bj) = 1. _d 0
214            ELSE            ELSE
215             pC(I,J,bi,bj) =  1. _d 0 / aC             pC(I,J,bi,bj) =  1. _d 0 / aC
216            ENDIF            ENDIF
# Line 177  C--   Update overlap regions Line 238  C--   Update overlap regions
238        _EXCH_XY_R4(pW, myThid)        _EXCH_XY_R4(pW, myThid)
239        _EXCH_XY_R4(pS, myThid)        _EXCH_XY_R4(pS, myThid)
240  CcnhDebugStarts  CcnhDebugStarts
241        CALL PLOT_FIELD_XYRS( pC, 'pC   INI_CG2D.2' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( pC, 'pC   INI_CG2D.2' , 1, myThid )
242        CALL PLOT_FIELD_XYRS( pW, 'pW   INI_CG2D.2' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( pW, 'pW   INI_CG2D.2' , 1, myThid )
243        CALL PLOT_FIELD_XYRS( pS, 'pS   INI_CG2D.2' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( pS, 'pS   INI_CG2D.2' , 1, myThid )
244  CcnhDebugEnds  CcnhDebugEnds
245    
 C--   Set default values for initial guess and RHS  
       IF ( startTime .EQ. 0 ) THEN  
        DO bj=myByLo(myThid),myByHi(myThid)  
         DO bi=myBxLo(myThid),myBxHi(myThid)  
          DO J=1,sNy  
           DO I=1,sNx  
            cg2d_x(I,J,bi,bj) = 0. _d 0  
            cg2d_b(I,J,bi,bj) = 0. _d 0  
 #ifdef ALLOW_CD  
            cg2d_xNM1(I,J,bi,bj) = 0. _d 0  
 #endif  
           ENDDO  
          ENDDO  
         ENDDO  
        ENDDO  
 C--    Update overlap regions  
        _EXCH_XY_R8(cg2d_x, myThid)  
        _EXCH_XY_R8(cg2d_b, myThid)  
 #ifdef ALLOW_CD  
        _EXCH_XY_R8(cg2d_xNM1, myThid)  
 #endif  
       ENDIF  
   
246        RETURN        RETURN
247        END        END
Legend:
Removed from v.1.18  
changed lines
  Added in v.1.31
  ViewVC Help
Powered by ViewVC 1.1.22