/[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.6 by cnh, Mon May 25 16:58:49 1998 UTC revision 1.35 by adcroft, Thu Jun 7 15:53:21 2001 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_EEOPTIONS.h"  #include "CPP_OPTIONS.h"
5    
 CStartOfInterface  
6        SUBROUTINE INI_CG2D( myThid )        SUBROUTINE INI_CG2D( myThid )
7  C     /==========================================================\  C     /==========================================================\
8  C     | SUBROUTINE INI_CG2D                                      |  C     | SUBROUTINE INI_CG2D                                      |
# Line 11  C     |================================= Line 11  C     |=================================
11  C     | These arrays are purely a function of the basin geom.    |  C     | These arrays are purely a function of the basin geom.    |
12  C     | We set then here once and them use then repeatedly.      |  C     | We set then here once and them use then repeatedly.      |
13  C     \==========================================================/  C     \==========================================================/
14          IMPLICIT NONE
15    
16  C     === Global variables ===  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"  #include "GRID.h"
21    c#include "DYNVARS.h"
22    #include "SURFACE.h"
23  #include "CG2D.h"  #include "CG2D.h"
24    #ifdef ALLOW_OBCS
25    #include "OBCS.h"
26    #endif
27    
28  C     === Routine arguments ===  C     === Routine arguments ===
29  C     myThid - Thread no. that called this routine.  C     myThid - Thread no. that called this routine.
30        INTEGER myThid        INTEGER myThid
 CEndOfInterface  
31    
32  C     === Local variables ===  C     === Local variables ===
33  C     xG, yG - Global coordinate location.  C     xG, yG - Global coordinate location.
# Line 31  C     iG, jG - Global coordinate index Line 36  C     iG, jG - Global coordinate index
36  C     bi,bj  - Loop counters  C     bi,bj  - Loop counters
37  C     faceArea - Temporary used to hold cell face areas.  C     faceArea - Temporary used to hold cell face areas.
38  C     I,J,K  C     I,J,K
39  C     myNorm - Work variable used in clculating normalisation factor  C     myNorm - Work variable used in calculating normalisation factor
40    C     sumArea - Work variable used to compute the total Domain Area
41        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
42        INTEGER bi, bj        INTEGER bi, bj
43        INTEGER I,  J, K        INTEGER I,  J, K
44        real    faceArea        _RL     faceArea
45        _RL     myNorm        _RS     myNorm
46          _RL     sumArea
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    C-  gBaro = gravity (except for External mode test with reduced gravity)
69            DO bj=myByLo(myThid),myByHi(myThid)
70             DO bi=myBxLo(myThid),myBxHi(myThid)
71              DO J=1-Oly,sNy+Oly
72               DO I=1-Olx,sNx+Olx
73                 Bo_surf(I,J,bi,bj) = gBaro
74                 recip_Bo(I,J,bi,bj) = 1. _d 0 / gBaro
75               ENDDO
76              ENDDO
77             ENDDO
78            ENDDO
79          ENDIF
80    
81    C--   Update overlap regions
82          _EXCH_XY_R8(Bo_surf, myThid)
83          _EXCH_XY_R8(recip_Bo, myThid)
84    
85  C--   Initialise laplace operator  C--   Initialise laplace operator
86  C     aW2d: integral in Z Ax/dX  C     aW2d: integral in Z Ax/dX
87  C     aS2d: integral in Z Ay/dY  C     aS2d: integral in Z Ay/dY
# Line 51  C     aS2d: integral in Z Ay/dY Line 94  C     aS2d: integral in Z Ay/dY
94            aS2d(I,J,bi,bj) = 0. _d 0            aS2d(I,J,bi,bj) = 0. _d 0
95           ENDDO           ENDDO
96          ENDDO          ENDDO
97          DO K=1,Nz          DO K=1,Nr
98           DO J=1,sNy           DO J=1,sNy
99            DO I=1,sNx            DO I=1,sNx
100             faceArea = dyG(I,J,bi,bj)*dzF(K)*HFacW(I,J,K,bi,bj)             faceArea = _dyG(I,J,bi,bj)*drF(K)
101         &               *_hFacW(I,J,K,bi,bj)
102             aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) +             aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) +
103       &      gBaro*faceArea*rDxC(I,J,bi,bj)       &      implicSurfPress*implicDiv2DFlow
104             faceArea = dxG(I,J,bi,bj)*dzF(K)*HFacS(I,J,K,bi,bj)       &           *faceArea*recip_dxC(I,J,bi,bj)
105               faceArea = _dxG(I,J,bi,bj)*drF(K)
106         &               *_hFacS(I,J,K,bi,bj)
107             aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) +             aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) +
108       &      gBaro*faceArea*rDyC(I,J,bi,bj)       &      implicSurfPress*implicDiv2DFlow
109         &           *faceArea*recip_dyC(I,J,bi,bj)
110            ENDDO            ENDDO
111           ENDDO           ENDDO
112          ENDDO          ENDDO
113    #ifdef ALLOW_OBCS
114            IF (useOBCS) THEN
115             DO I=1,sNx
116              IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj),bi,bj)=0.
117              IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj)+1,bi,bj)=0.
118              IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj)+1,bi,bj)=0.
119              IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj),bi,bj)=0.
120             ENDDO
121             DO J=1,sNy
122              IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj),J,bi,bj)=0.
123              IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj)+1,J,bi,bj)=0.
124              IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj)+1,J,bi,bj)=0.
125              IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj),J,bi,bj)=0.
126             ENDDO
127            ENDIF
128    #endif
129          DO J=1,sNy          DO J=1,sNy
130           DO I=1,sNx           DO I=1,sNx
131            myNorm = MAX(ABS(aW2d(I,J,bi,bj)),myNorm)            myNorm = MAX(ABS(aW2d(I,J,bi,bj)),myNorm)
# Line 71  C     aS2d: integral in Z Ay/dY Line 134  C     aS2d: integral in Z Ay/dY
134          ENDDO          ENDDO
135         ENDDO         ENDDO
136        ENDDO        ENDDO
137        cg2dNbuf(1,myThid) = myNorm        _GLOBAL_MAX_R4( myNorm, myThid )
138        _GLOBAL_MAX_R8( cg2dNbuf, myNorm, myThid )        IF ( myNorm .NE. 0. _d 0 ) THEN
139        IF ( cg2dNbuf(1,1) .NE. 0. _d 0 ) THEN         myNorm = 1. _d 0/myNorm
        myNorm = 1. _d 0/cg2dNbuf(1,1)  
140        ELSE        ELSE
141         myNorm = 1. _d 0         myNorm = 1. _d 0
142        ENDIF        ENDIF
143         cg2dNorm = myNorm         cg2dNorm = myNorm
144        _BEGIN_MASTER( myThid )        _BEGIN_MASTER( myThid )
145  CcnhDebugStarts  CcnhDebugStarts
146         WRITE(msgBuf,'(A,F)') '// CG2D normalisation factor = ', cg2dNorm         WRITE(msgBuf,'(A,E40.25)') '// CG2D normalisation factor = ',
147         &               cg2dNorm
148         CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)         CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
149         WRITE(msgBuf,*) '                               '         WRITE(msgBuf,*) '                               '
150         CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)         CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
# Line 100  CcnhDebugEnds Line 163  CcnhDebugEnds
163                
164  C--   Update overlap regions  C--   Update overlap regions
165  CcnhDebugStarts  CcnhDebugStarts
166  C     CALL PLOT_FIELD_XYR8( aW2d, 'AW2D INI_CG2D.1' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.1' , 1, myThid )
167  C     CALL PLOT_FIELD_XYR8( aS2d, 'AS2D INI_CG2D.1' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.1' , 1, myThid )
168  CcnhDebugEnds  CcnhDebugEnds
169        _EXCH_XY_R4(aW2d, myThid)  c     _EXCH_XY_R4(aW2d, myThid)
170        _EXCH_XY_R4(aS2d, myThid)  c     _EXCH_XY_R4(aS2d, myThid)
171          CALL EXCH_UV_XY_RS(aW2d,aS2d,.FALSE.,myThid)
172  CcnhDebugStarts  CcnhDebugStarts
173  C     CALL PLOT_FIELD_XYR8( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )
174  C     CALL PLOT_FIELD_XYR8( aS2d, 'AS2D INI_CG2D.2' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.2' , 1, myThid )
175  CcnhDebugEnds  CcnhDebugEnds
176    
177    C-- Define the solver tolerance in the appropriate Unit :
178          cg2dNormaliseRHS = cg2dTargetResWunit.LE.0
179          IF (cg2dNormaliseRHS) THEN
180    C-  when using a normalisation of RHS, tolerance has no unit => no conversion
181            cg2dTolerance = cg2dTargetResidual
182          ELSE
183    C-  compute the total Area of the domain :
184            sumArea = 0.
185            DO bj=myByLo(myThid),myByHi(myThid)
186             DO bi=myBxLo(myThid),myBxHi(myThid)
187              DO j=1,sNy
188               DO i=1,sNx
189                 IF (Ro_surf(i,j,bi,bj).GT.R_low(i,j,bi,bj)) THEN
190                   sumArea = sumArea + rA(i,j,bi,bj)
191                 ENDIF
192               ENDDO
193              ENDDO
194             ENDDO
195            ENDDO
196    c       WRITE(*,*) ' mythid, sumArea = ', mythid, sumArea
197            _GLOBAL_SUM_R8( sumArea, myThid )
198    C-  convert Target-Residual (in W unit) to cg2d-solver residual unit [m^2/s^2]
199            cg2dTolerance = cg2dNorm * cg2dTargetResWunit
200         &                           * sumArea / deltaTMom
201            WRITE(*,'(2A,1P2E22.14)') ' ini_cg2d: ',
202         &          'sumArea,cg2dTolerance =', sumArea,cg2dTolerance
203          ENDIF
204        
205  C--   Initialise preconditioner  C--   Initialise preconditioner
206  C     Note. 20th May 1998  C     Note. 20th May 1998
207  C           I made a weird discovery! In the model paper we argue  C           I made a weird discovery! In the model paper we argue
# Line 131  C           defaults to 0.51 but can be Line 223  C           defaults to 0.51 but can be
223            pC(I,J,bi,bj) = 1. _d 0            pC(I,J,bi,bj) = 1. _d 0
224            aC = -(            aC = -(
225       &     aW2d(I,J,bi,bj) + aW2d(I+1,J  ,bi,bj)       &     aW2d(I,J,bi,bj) + aW2d(I+1,J  ,bi,bj)
226       &    +aS2d(I,J,bi,bj) + aS2D(I  ,J+1,bi,bj)       &    +aS2d(I,J,bi,bj) + aS2d(I  ,J+1,bi,bj)
227       &    +freeSurfFac*myNorm*       &    +freeSurfFac*myNorm*recip_Bo(I,J,bi,bj)*
228       &     zA(I,J,bi,bj)/deltaTMom/deltaTMom       &     rA(I,J,bi,bj)/deltaTMom/deltaTMom
229       &    )       &    )
230            aCs = -(            aCs = -(
231       &     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)
232       &    +aS2d(I,J-1,bi,bj) + aS2d(I  ,J  ,bi,bj)       &    +aS2d(I,J-1,bi,bj) + aS2d(I  ,J  ,bi,bj)
233       &    +freeSurfFac*myNorm*       &    +freeSurfFac*myNorm*recip_Bo(I,J-1,bi,bj)*
234       &     zA(I,J-1,bi,bj)/deltaTMom/deltaTMom       &     rA(I,J-1,bi,bj)/deltaTMom/deltaTMom
235       &    )       &    )
236            aCw = -(            aCw = -(
237       &     aW2d(I-1,J,bi,bj) + aW2d(I  ,J  ,bi,bj)       &     aW2d(I-1,J,bi,bj) + aW2d(I  ,J  ,bi,bj)
238       &    +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)
239       &    +freeSurfFac*myNorm*       &    +freeSurfFac*myNorm*recip_Bo(I-1,J,bi,bj)*
240       &     zA(I-1,J,bi,bj)/deltaTMom/deltaTMom       &     rA(I-1,J,bi,bj)/deltaTMom/deltaTMom
241       &    )       &    )
242            IF ( aC .EQ. 0. ) THEN            IF ( aC .EQ. 0. ) THEN
243              pC(I,J,bi,bj) = 0. _d 0              pC(I,J,bi,bj) = 1. _d 0
244            ELSE            ELSE
245             pC(I,J,bi,bj) =  1. _d 0 / aC             pC(I,J,bi,bj) =  1. _d 0 / aC
246            ENDIF            ENDIF
# Line 173  C         pS(I,J,bi,bj) = 0. Line 265  C         pS(I,J,bi,bj) = 0.
265        ENDDO        ENDDO
266  C--   Update overlap regions  C--   Update overlap regions
267        _EXCH_XY_R4(pC, myThid)        _EXCH_XY_R4(pC, myThid)
268        _EXCH_XY_R4(pW, myThid)  c     _EXCH_XY_R4(pW, myThid)
269        _EXCH_XY_R4(pS, myThid)  c     _EXCH_XY_R4(pS, myThid)
270          CALL EXCH_UV_XY_RS(pW,pS,.FALSE.,myThid)
271  C--   Set default values for initial guess and RHS  CcnhDebugStarts
272        IF ( startTime .EQ. 0 ) THEN  C     CALL PLOT_FIELD_XYRS( pC, 'pC   INI_CG2D.2' , 1, myThid )
273         DO bj=myByLo(myThid),myByHi(myThid)  C     CALL PLOT_FIELD_XYRS( pW, 'pW   INI_CG2D.2' , 1, myThid )
274          DO bi=myBxLo(myThid),myBxHi(myThid)  C     CALL PLOT_FIELD_XYRS( pS, 'pS   INI_CG2D.2' , 1, myThid )
275           DO J=1,sNy  CcnhDebugEnds
           DO I=1,sNx  
            cg2d_x(I,J,bi,bj) = 0. _d 0  
            cg2d_b(I,J,bi,bj) = 0. _d 0  
           ENDDO  
          ENDDO  
         ENDDO  
        ENDDO  
 C--    Update overlap regions  
        _EXCH_XY_R8(cg2d_x, myThid)  
        _EXCH_XY_R8(cg2d_b, myThid)  
       ENDIF  
276    
277        RETURN        RETURN
278        END        END
Legend:
Removed from v.1.6  
changed lines
  Added in v.1.35
  ViewVC Help
Powered by ViewVC 1.1.22