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

Annotation of /MITgcm/model/src/update_cg2d.F

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

Revision 1.2 - (hide annotations) (download)
Wed Sep 26 18:09:16 2001 UTC (22 years, 8 months ago) by cnh
Branch: MAIN
CVS Tags: checkpoint43a-release1mods, release1_b1, checkpoint43, icebear5, icebear4, icebear3, icebear2, release1-branch_tutorials, chkpt44a_post, chkpt44c_pre, ecco_c44_e19, ecco_c44_e18, ecco_c44_e17, ecco_c44_e16, release1-branch-end, checkpoint44b_post, ecco_ice2, ecco_ice1, ecco_c44_e22, ecco_c44_e25, chkpt44a_pre, ecco_c44_e23, ecco_c44_e20, ecco_c44_e21, ecco_c44_e26, ecco_c44_e27, ecco_c44_e24, ecco-branch-mod1, ecco-branch-mod2, ecco-branch-mod3, ecco-branch-mod4, ecco-branch-mod5, release1_beta1, checkpoint44b_pre, checkpoint42, checkpoint41, checkpoint44, chkpt44c_post, release1-branch_branchpoint
Branch point for: c24_e25_ice, release1-branch, release1, ecco-branch, icebear, release1_coupled
Changes since 1.1: +21 -12 lines 
Bringing comments up to data and formatting for document extraction.
1 cnh 1.2 C $Header: /u/gcmpack/models/MITgcmUV/model/src/update_cg2d.F,v 1.1 2001/08/27 18:50:10 jmc Exp $
2 jmc 1.1 C $Name: $
3    
4     #include "CPP_OPTIONS.h"
5 cnh 1.2 CBOP
6     C !ROUTINE: UPDATE_CG2D
7     C !INTERFACE:
8     SUBROUTINE UPDATE_CG2D( myTime, myIter, myThid )
9     C !DESCRIPTION: \bv
10     C *==========================================================*
11     C | SUBROUTINE UPDATE_CG2D
12     C | o Update 2d conjugate gradient solver operators
13     C | account for Free-Surf effect on total column thickness
14     C *==========================================================*
15     C | This routine is based on INI_CG2D, and simplified. It is
16     C | only active when the non-linear free surface mode of
17     C | equations is active.
18     C *==========================================================*
19     C \ev
20 jmc 1.1
21 cnh 1.2 C !USES:
22 jmc 1.1 IMPLICIT NONE
23     C === Global variables ===
24     #include "SIZE.h"
25     #include "EEPARAMS.h"
26     #include "PARAMS.h"
27     #include "GRID.h"
28     #include "SURFACE.h"
29     #include "CG2D.h"
30     #ifdef ALLOW_OBCS
31     #include "OBCS.h"
32     #endif
33    
34 cnh 1.2 C !INPUT/OUTPUT PARAMETERS:
35 jmc 1.1 C === Routine arguments ===
36     C myTime - Current time in simulation
37     C myIter - Current iteration number in simulation
38     C myThid - Thread number for this instance of the routine.
39     _RL myTime
40     INTEGER myIter
41     INTEGER myThid
42 cnh 1.2
43     C !LOCAL VARIABLES:
44 jmc 1.1 #ifdef NONLIN_FRSURF
45     C-- Note : compared to "INI_CG2D", no needs to compute again
46     C the solver norn=malisation factor of the solver tolerance
47     C === Local variables ===
48     C bi,bj - Loop counters
49     C I,J,K
50     C faceArea - Temporary used to hold cell face areas.
51     INTEGER bi, bj
52     INTEGER I, J, K
53     _RL faceArea
54     _RL aC, aCw, aCs
55 cnh 1.2 CEOP
56 jmc 1.1
57     C-- Initialise laplace operator
58     C aW2d: integral in Z Ax/dX
59     C aS2d: integral in Z Ay/dY
60     DO bj=myByLo(myThid),myByHi(myThid)
61     DO bi=myBxLo(myThid),myBxHi(myThid)
62     DO J=1,sNy
63     DO I=1,sNx
64     aW2d(I,J,bi,bj) = 0. _d 0
65     aS2d(I,J,bi,bj) = 0. _d 0
66     ENDDO
67     ENDDO
68     DO K=1,Nr
69     DO J=1,sNy
70     DO I=1,sNx
71     faceArea = _dyG(I,J,bi,bj)*drF(K)
72     & *_hFacW(I,J,K,bi,bj)
73     aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) +
74     & faceArea*recip_dxC(I,J,bi,bj)
75     faceArea = _dxG(I,J,bi,bj)*drF(K)
76     & *_hFacS(I,J,K,bi,bj)
77     aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) +
78     & faceArea*recip_dyC(I,J,bi,bj)
79     ENDDO
80     ENDDO
81     ENDDO
82     #ifdef ALLOW_OBCS
83     IF (useOBCS) THEN
84     DO I=1,sNx
85     IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj),bi,bj)=0.
86     IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj)+1,bi,bj)=0.
87     IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj)+1,bi,bj)=0.
88     IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj),bi,bj)=0.
89     ENDDO
90     DO J=1,sNy
91     IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj),J,bi,bj)=0.
92     IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj)+1,J,bi,bj)=0.
93     IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj)+1,J,bi,bj)=0.
94     IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj),J,bi,bj)=0.
95     ENDDO
96     ENDIF
97     #endif
98     DO J=1,sNy
99     DO I=1,sNx
100     aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj)*cg2dNorm
101     & *implicSurfPress*implicDiv2DFlow
102     aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj)*cg2dNorm
103     & *implicSurfPress*implicDiv2DFlow
104     ENDDO
105     ENDDO
106     ENDDO
107     ENDDO
108    
109     C-- Update overlap regions
110     c _EXCH_XY_R4(aW2d, myThid)
111     c _EXCH_XY_R4(aS2d, myThid)
112     CALL EXCH_UV_XY_RS(aW2d,aS2d,.FALSE.,myThid)
113    
114     C-- Initialise preconditioner
115     DO bj=myByLo(myThid),myByHi(myThid)
116     DO bi=myBxLo(myThid),myBxHi(myThid)
117     DO J=1,sNy
118     DO I=1,sNx
119     pC(I,J,bi,bj) = 1. _d 0
120     aC = -(
121     & aW2d(I,J,bi,bj) + aW2d(I+1,J ,bi,bj)
122     & +aS2d(I,J,bi,bj) + aS2d(I ,J+1,bi,bj)
123     & +freeSurfFac*cg2dNorm*recip_Bo(I,J,bi,bj)*
124     & rA(I,J,bi,bj)/deltaTMom/deltaTMom
125     & )
126     aCs = -(
127     & aW2d(I,J-1,bi,bj) + aW2d(I+1,J-1,bi,bj)
128     & +aS2d(I,J-1,bi,bj) + aS2d(I ,J ,bi,bj)
129     & +freeSurfFac*cg2dNorm*recip_Bo(I,J-1,bi,bj)*
130     & rA(I,J-1,bi,bj)/deltaTMom/deltaTMom
131     & )
132     aCw = -(
133     & aW2d(I-1,J,bi,bj) + aW2d(I ,J ,bi,bj)
134     & +aS2d(I-1,J,bi,bj) + aS2d(I-1,J+1,bi,bj)
135     & +freeSurfFac*cg2dNorm*recip_Bo(I-1,J,bi,bj)*
136     & rA(I-1,J,bi,bj)/deltaTMom/deltaTMom
137     & )
138     IF ( aC .EQ. 0. ) THEN
139     pC(I,J,bi,bj) = 1. _d 0
140     ELSE
141     pC(I,J,bi,bj) = 1. _d 0 / aC
142     ENDIF
143     IF ( aC + aCw .EQ. 0. ) THEN
144     pW(I,J,bi,bj) = 0.
145     ELSE
146     pW(I,J,bi,bj) =
147     & -aW2d(I ,J ,bi,bj)/((cg2dpcOffDFac *(aCw+aC))**2 )
148     ENDIF
149     IF ( aC + aCs .EQ. 0. ) THEN
150     pS(I,J,bi,bj) = 0.
151     ELSE
152     pS(I,J,bi,bj) =
153     & -aS2d(I ,J ,bi,bj)/((cg2dpcOffDFac *(aCs+aC))**2 )
154     ENDIF
155     ENDDO
156     ENDDO
157     ENDDO
158     ENDDO
159     C-- Update overlap regions
160     _EXCH_XY_R4(pC, myThid)
161     c _EXCH_XY_R4(pW, myThid)
162     c _EXCH_XY_R4(pS, myThid)
163     CALL EXCH_UV_XY_RS(pW,pS,.FALSE.,myThid)
164    
165     #endif /* NONLIN_FRSURF */
166    
167     RETURN
168     END
  ViewVC Help
Powered by ViewVC 1.1.22