/[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.38 by jmc, Sun Feb 10 20:04:11 2002 UTC revision 1.44 by jmc, Tue Oct 17 18:52:34 2006 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4    #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    
7  CBOP  CBOP
# Line 25  C     === Global variables === Line 26  C     === Global variables ===
26  #include "EEPARAMS.h"  #include "EEPARAMS.h"
27  #include "PARAMS.h"  #include "PARAMS.h"
28  #include "GRID.h"  #include "GRID.h"
 c#include "DYNVARS.h"  
29  #include "SURFACE.h"  #include "SURFACE.h"
30  #include "CG2D.h"  #include "CG2D.h"
31  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
# Line 52  C     sumArea - Work variable used to co Line 52  C     sumArea - Work variable used to co
52        INTEGER I,  J, K        INTEGER I,  J, K
53        _RL     faceArea        _RL     faceArea
54        _RS     myNorm        _RS     myNorm
       _RL     sumArea  
55        _RL     aC, aCw, aCs        _RL     aC, aCw, aCs
56  CEOP  CEOP
57    
# Line 137  C     aS2d: integral in Z Ay/dY Line 136  C     aS2d: integral in Z Ay/dY
136        ELSE        ELSE
137         myNorm = 1. _d 0         myNorm = 1. _d 0
138        ENDIF        ENDIF
        cg2dNorm = myNorm  
       _BEGIN_MASTER( myThid )  
 CcnhDebugStarts  
        WRITE(msgBuf,'(A,E40.25)') '// CG2D normalisation factor = ',  
      &               cg2dNorm  
        CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)  
        WRITE(msgBuf,*) '                               '  
        CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)  
 CcnhDebugEnds  
       _END_MASTER( myThid )  
139        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
140         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
141          DO J=1,sNy          DO J=1,sNy
# Line 171  C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D Line 160  C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D
160  C     CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.2' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.2' , 1, myThid )
161  CcnhDebugEnds  CcnhDebugEnds
162    
163          _BEGIN_MASTER(myThid)
164    C-- set global parameter in common block:
165          cg2dNorm = myNorm
166  C-- Define the solver tolerance in the appropriate Unit :  C-- Define the solver tolerance in the appropriate Unit :
167        cg2dNormaliseRHS = cg2dTargetResWunit.LE.0        cg2dNormaliseRHS = cg2dTargetResWunit.LE.0.
168        IF (cg2dNormaliseRHS) THEN        IF (cg2dNormaliseRHS) THEN
169  C-  when using a normalisation of RHS, tolerance has no unit => no conversion  C-  when using a normalisation of RHS, tolerance has no unit => no conversion
170          cg2dTolerance = cg2dTargetResidual          cg2dTolerance = cg2dTargetResidual
171        ELSE        ELSE
 C-  compute the total Area of the domain :  
         sumArea = 0.  
         DO bj=myByLo(myThid),myByHi(myThid)  
          DO bi=myBxLo(myThid),myBxHi(myThid)  
           DO j=1,sNy  
            DO i=1,sNx  
              IF (Ro_surf(i,j,bi,bj).GT.R_low(i,j,bi,bj)) THEN  
                sumArea = sumArea + rA(i,j,bi,bj)  
              ENDIF  
            ENDDO  
           ENDDO  
          ENDDO  
         ENDDO  
 c       WRITE(*,*) ' mythid, sumArea = ', mythid, sumArea  
         _GLOBAL_SUM_R8( sumArea, myThid )  
172  C-  convert Target-Residual (in W unit) to cg2d-solver residual unit [m^2/s^2]  C-  convert Target-Residual (in W unit) to cg2d-solver residual unit [m^2/s^2]
173          cg2dTolerance = cg2dNorm * cg2dTargetResWunit          cg2dTolerance = cg2dNorm * cg2dTargetResWunit
174       &                           * sumArea / deltaTMom       &                           * globalArea / deltaTmom
         WRITE(*,'(2A,1P2E22.14)') ' ini_cg2d: ',  
      &          'sumArea,cg2dTolerance =', sumArea,cg2dTolerance  
175        ENDIF        ENDIF
176          _END_MASTER(myThid)
177    
178    CcnhDebugStarts
179          _BEGIN_MASTER( myThid )
180           WRITE(msgBuf,'(2A,1PE23.16)') 'INI_CG2D: ',
181         &      'CG2D normalisation factor = ', cg2dNorm
182           CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
183           IF (.NOT.cg2dNormaliseRHS) THEN
184            WRITE(msgBuf,'(2A,1PE22.15,A,1PE16.10,A)') 'INI_CG2D: ',
185         &      'cg2dTolerance =', cg2dTolerance, ' (Area=',globalArea,')'
186            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
187           ENDIF
188           WRITE(msgBuf,*) ' '
189           CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
190          _END_MASTER( myThid )
191    CcnhDebugEnds
192            
193  C--   Initialise preconditioner  C--   Initialise preconditioner
194  C     Note. 20th May 1998  C     Note. 20th May 1998
# Line 222  C           defaults to 0.51 but can be Line 213  C           defaults to 0.51 but can be
213       &     aW2d(I,J,bi,bj) + aW2d(I+1,J  ,bi,bj)       &     aW2d(I,J,bi,bj) + aW2d(I+1,J  ,bi,bj)
214       &    +aS2d(I,J,bi,bj) + aS2d(I  ,J+1,bi,bj)       &    +aS2d(I,J,bi,bj) + aS2d(I  ,J+1,bi,bj)
215       &    +freeSurfFac*myNorm*recip_Bo(I,J,bi,bj)*       &    +freeSurfFac*myNorm*recip_Bo(I,J,bi,bj)*
216       &     rA(I,J,bi,bj)/deltaTMom/deltaTMom       &     rA(I,J,bi,bj)/deltaTMom/deltaTfreesurf
217       &    )       &    )
218            aCs = -(            aCs = -(
219       &     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)
220       &    +aS2d(I,J-1,bi,bj) + aS2d(I  ,J  ,bi,bj)       &    +aS2d(I,J-1,bi,bj) + aS2d(I  ,J  ,bi,bj)
221       &    +freeSurfFac*myNorm*recip_Bo(I,J-1,bi,bj)*       &    +freeSurfFac*myNorm*recip_Bo(I,J-1,bi,bj)*
222       &     rA(I,J-1,bi,bj)/deltaTMom/deltaTMom       &     rA(I,J-1,bi,bj)/deltaTMom/deltaTfreesurf
223       &    )       &    )
224            aCw = -(            aCw = -(
225       &     aW2d(I-1,J,bi,bj) + aW2d(I  ,J  ,bi,bj)       &     aW2d(I-1,J,bi,bj) + aW2d(I  ,J  ,bi,bj)
226       &    +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)
227       &    +freeSurfFac*myNorm*recip_Bo(I-1,J,bi,bj)*       &    +freeSurfFac*myNorm*recip_Bo(I-1,J,bi,bj)*
228       &     rA(I-1,J,bi,bj)/deltaTMom/deltaTMom       &     rA(I-1,J,bi,bj)/deltaTMom/deltaTfreesurf
229       &    )       &    )
230            IF ( aC .EQ. 0. ) THEN            IF ( aC .EQ. 0. ) THEN
231              pC(I,J,bi,bj) = 1. _d 0              pC(I,J,bi,bj) = 1. _d 0
Legend:
Removed from v.1.38  
changed lines
  Added in v.1.44
  ViewVC Help
Powered by ViewVC 1.1.22