/[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.33.2.3 by adcroft, Tue Apr 3 18:00:50 2001 UTC revision 1.41 by jmc, Mon Jan 13 01:31:22 2003 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    
6    CBOP
7    C     !ROUTINE: INI_CG2D
8    C     !INTERFACE:
9        SUBROUTINE INI_CG2D( myThid )        SUBROUTINE INI_CG2D( myThid )
 C     /==========================================================\  
 C     | SUBROUTINE INI_CG2D                                      |  
 C     | o Initialise 2d conjugate gradient solver operators.     |  
 C     |==========================================================|  
 C     | These arrays are purely a function of the basin geom.    |  
 C     | We set then here once and them use then repeatedly.      |  
 C     \==========================================================/  
       IMPLICIT NONE  
10    
11    C     !DESCRIPTION: \bv
12    C     *==========================================================*
13    C     | SUBROUTINE INI_CG2D                                      
14    C     | o Initialise 2d conjugate gradient solver operators.      
15    C     *==========================================================*
16    C     | These arrays are purely a function of the basin geom.    
17    C     | We set then here once and them use then repeatedly.      
18    C     *==========================================================*
19    C     \ev
20    
21    C     !USES:
22          IMPLICIT NONE
23  C     === Global variables ===  C     === Global variables ===
24  #include "SIZE.h"  #include "SIZE.h"
25  #include "EEPARAMS.h"  #include "EEPARAMS.h"
# Line 25  c#include "DYNVARS.h" Line 32  c#include "DYNVARS.h"
32  #include "OBCS.h"  #include "OBCS.h"
33  #endif  #endif
34    
35    C     !INPUT/OUTPUT PARAMETERS:
36  C     === Routine arguments ===  C     === Routine arguments ===
37  C     myThid - Thread no. that called this routine.  C     myThid - Thread no. that called this routine.
38        INTEGER myThid        INTEGER myThid
39    
40    C     !LOCAL VARIABLES:
41  C     === Local variables ===  C     === Local variables ===
42  C     xG, yG - Global coordinate location.  C     xG, yG - Global coordinate location.
43  C     zG  C     zG
# Line 36  C     iG, jG - Global coordinate index Line 45  C     iG, jG - Global coordinate index
45  C     bi,bj  - Loop counters  C     bi,bj  - Loop counters
46  C     faceArea - Temporary used to hold cell face areas.  C     faceArea - Temporary used to hold cell face areas.
47  C     I,J,K  C     I,J,K
48  C     myNorm - Work variable used in clculating normalisation factor  C     myNorm - Work variable used in calculating normalisation factor
49    C     sumArea - Work variable used to compute the total Domain Area
50        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
51        INTEGER bi, bj        INTEGER bi, bj
52        INTEGER I,  J, K        INTEGER I,  J, K
53        _RL     faceArea        _RL     faceArea
54        _RS     myNorm        _RS     myNorm
55          _RL     sumArea
56        _RL     aC, aCw, aCs        _RL     aC, aCw, aCs
57    CEOP
58    
59  C-- Initialise -Boyancy at surface level : Bo_surf  C--   Initialize arrays in common blocs (CG2D.h) ; not really necessary
60  C   Bo_surf is defined as d/dr(Phi_surf) and set to g/rtoz (linear free surface)  C     but safer when EXCH do not fill all the overlap regions.
61  C     with rtoz = conversion factor from r-unit to z-unit  (=horiVertRatio)        DO bj=myByLo(myThid),myByHi(myThid)
62  C   an acurate formulation will include P_surf and T,S_surf effects on rho_surf:         DO bi=myBxLo(myThid),myBxHi(myThid)
63  C    z-ocean (rtoz=1) : Bo_surf = - Boyancy = gravity * rho_surf/rho_0          DO J=1-OLy,sNy+OLy
64  C    p-atmos (rtoz=rho_c*g) : Bo_surf = (1/rho)_surf           DO I=1-OLx,sNx+OLx
65  C--            aW2d(I,J,bi,bj) = 0. _d 0
66        IF ( buoyancyRelation .EQ. 'ATMOSPHERIC' ) THEN            aS2d(I,J,bi,bj) = 0. _d 0
67          DO bj=myByLo(myThid),myByHi(myThid)            pW(I,J,bi,bj) = 0. _d 0
68           DO bi=myBxLo(myThid),myBxHi(myThid)            pS(I,J,bi,bj) = 0. _d 0
69            DO J=1-Oly,sNy+Oly            pC(I,J,bi,bj) = 0. _d 0
70             DO I=1-Olx,sNx+Olx            cg2d_q(I,J,bi,bj) = 0. _d 0
              Bo_surf(I,J,bi,bj) = recip_rhoConst  
              recip_Bo(I,J,bi,bj) = rhoConst  
            ENDDO  
           ENDDO  
71           ENDDO           ENDDO
72          ENDDO          ENDDO
73        ELSE          DO J=1-1,sNy+1
74  C-  gBaro = gravity (except for External mode test with reduced gravity)           DO I=1-1,sNx+1
75          DO bj=myByLo(myThid),myByHi(myThid)            cg2d_r(I,J,bi,bj) = 0. _d 0
76           DO bi=myBxLo(myThid),myBxHi(myThid)            cg2d_s(I,J,bi,bj) = 0. _d 0
           DO J=1-Oly,sNy+Oly  
            DO I=1-Olx,sNx+Olx  
              Bo_surf(I,J,bi,bj) = gBaro  
              recip_Bo(I,J,bi,bj) = 1. _d 0 / gBaro  
            ENDDO  
           ENDDO  
77           ENDDO           ENDDO
78          ENDDO          ENDDO
79        ENDIF         ENDDO
80          ENDDO
 C--   Update overlap regions  
       _EXCH_XY_R8(Bo_surf, myThid)  
       _EXCH_XY_R8(recip_Bo, myThid)  
81    
82  C--   Initialise laplace operator  C--   Initialise laplace operator
83  C     aW2d: integral in Z Ax/dX  C     aW2d: integral in Z Ax/dX
# Line 172  C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D Line 171  C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D
171  C     CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.2' , 1, myThid )  C     CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.2' , 1, myThid )
172  CcnhDebugEnds  CcnhDebugEnds
173    
174    C-- Define the solver tolerance in the appropriate Unit :
175          cg2dNormaliseRHS = cg2dTargetResWunit.LE.0
176          IF (cg2dNormaliseRHS) THEN
177    C-  when using a normalisation of RHS, tolerance has no unit => no conversion
178            cg2dTolerance = cg2dTargetResidual
179          ELSE
180    C-  compute the total Area of the domain :
181            sumArea = 0.
182            DO bj=myByLo(myThid),myByHi(myThid)
183             DO bi=myBxLo(myThid),myBxHi(myThid)
184              DO j=1,sNy
185               DO i=1,sNx
186                 IF (Ro_surf(i,j,bi,bj).GT.R_low(i,j,bi,bj)) THEN
187                   sumArea = sumArea + rA(i,j,bi,bj)
188                 ENDIF
189               ENDDO
190              ENDDO
191             ENDDO
192            ENDDO
193    c       WRITE(*,*) ' mythid, sumArea = ', mythid, sumArea
194            _GLOBAL_SUM_R8( sumArea, myThid )
195    C-  convert Target-Residual (in W unit) to cg2d-solver residual unit [m^2/s^2]
196            cg2dTolerance = cg2dNorm * cg2dTargetResWunit
197         &                           * sumArea / deltaTmom
198            _BEGIN_MASTER( myThid )
199            WRITE(standardMessageUnit,'(2A,1P2E22.14)') ' ini_cg2d: ',
200         &          'sumArea,cg2dTolerance =', sumArea,cg2dTolerance
201            _END_MASTER( myThid )
202          ENDIF
203        
204  C--   Initialise preconditioner  C--   Initialise preconditioner
205  C     Note. 20th May 1998  C     Note. 20th May 1998
206  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 195  C           defaults to 0.51 but can be Line 224  C           defaults to 0.51 but can be
224       &     aW2d(I,J,bi,bj) + aW2d(I+1,J  ,bi,bj)       &     aW2d(I,J,bi,bj) + aW2d(I+1,J  ,bi,bj)
225       &    +aS2d(I,J,bi,bj) + aS2d(I  ,J+1,bi,bj)       &    +aS2d(I,J,bi,bj) + aS2d(I  ,J+1,bi,bj)
226       &    +freeSurfFac*myNorm*recip_Bo(I,J,bi,bj)*       &    +freeSurfFac*myNorm*recip_Bo(I,J,bi,bj)*
227       &     rA(I,J,bi,bj)/deltaTMom/deltaTMom       &     rA(I,J,bi,bj)/deltaTMom/deltaTfreesurf
228       &    )       &    )
229            aCs = -(            aCs = -(
230       &     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)
231       &    +aS2d(I,J-1,bi,bj) + aS2d(I  ,J  ,bi,bj)       &    +aS2d(I,J-1,bi,bj) + aS2d(I  ,J  ,bi,bj)
232       &    +freeSurfFac*myNorm*recip_Bo(I,J-1,bi,bj)*       &    +freeSurfFac*myNorm*recip_Bo(I,J-1,bi,bj)*
233       &     rA(I,J-1,bi,bj)/deltaTMom/deltaTMom       &     rA(I,J-1,bi,bj)/deltaTMom/deltaTfreesurf
234       &    )       &    )
235            aCw = -(            aCw = -(
236       &     aW2d(I-1,J,bi,bj) + aW2d(I  ,J  ,bi,bj)       &     aW2d(I-1,J,bi,bj) + aW2d(I  ,J  ,bi,bj)
237       &    +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)
238       &    +freeSurfFac*myNorm*recip_Bo(I-1,J,bi,bj)*       &    +freeSurfFac*myNorm*recip_Bo(I-1,J,bi,bj)*
239       &     rA(I-1,J,bi,bj)/deltaTMom/deltaTMom       &     rA(I-1,J,bi,bj)/deltaTMom/deltaTfreesurf
240       &    )       &    )
241            IF ( aC .EQ. 0. ) THEN            IF ( aC .EQ. 0. ) THEN
242              pC(I,J,bi,bj) = 1. _d 0              pC(I,J,bi,bj) = 1. _d 0
Legend:
Removed from v.1.33.2.3  
changed lines
  Added in v.1.41
  ViewVC Help
Powered by ViewVC 1.1.22