/[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.48 by mlosch, Mon Nov 23 16:15:54 2009 UTC revision 1.49 by jmc, Tue Dec 8 00:10:26 2009 UTC
# Line 40  C     myThid - Thread no. that called th Line 40  C     myThid - Thread no. that called th
40  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
41  C     === Local variables ===  C     === Local variables ===
42  C     bi,bj  :: tile indices  C     bi,bj  :: tile indices
43  C     I,J,K  :: Loop counters  C     i,j,k  :: Loop counters
44  C     faceArea - Temporary used to hold cell face areas.  C     faceArea - Temporary used to hold cell face areas.
45  C     myNorm - Work variable used in calculating normalisation factor  C     myNorm - Work variable used in calculating normalisation factor
46  C     sumArea - Work variable used to compute the total Domain Area  C     sumArea - Work variable used to compute the total Domain Area
47        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
48        INTEGER bi, bj        INTEGER bi, bj
49        INTEGER I, J, K, ks        INTEGER i, j, k, ks
50        _RL     faceArea        _RL     faceArea
51        _RS     myNorm        _RS     myNorm
52        _RS     aC, aCw, aCs        _RS     aC, aCw, aCs
# Line 56  C--   Initialize arrays in common blocs Line 56  C--   Initialize arrays in common blocs
56  C     but safer when EXCH do not fill all the overlap regions.  C     but safer when EXCH do not fill all the overlap regions.
57        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
58         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
59          DO J=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
60           DO I=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
61            aW2d(I,J,bi,bj) = 0. _d 0            aW2d(i,j,bi,bj) = 0. _d 0
62            aS2d(I,J,bi,bj) = 0. _d 0            aS2d(i,j,bi,bj) = 0. _d 0
63            aC2d(I,J,bi,bj) = 0. _d 0            aC2d(i,j,bi,bj) = 0. _d 0
64            pW(I,J,bi,bj) = 0. _d 0            pW(i,j,bi,bj) = 0. _d 0
65            pS(I,J,bi,bj) = 0. _d 0            pS(i,j,bi,bj) = 0. _d 0
66            pC(I,J,bi,bj) = 0. _d 0            pC(i,j,bi,bj) = 0. _d 0
67           ENDDO           ENDDO
68          ENDDO          ENDDO
69          DO J=1-1,sNy+1          DO j=1-1,sNy+1
70           DO I=1-1,sNx+1           DO i=1-1,sNx+1
71            cg2d_q(I,J,bi,bj) = 0. _d 0            cg2d_q(i,j,bi,bj) = 0. _d 0
72            cg2d_r(I,J,bi,bj) = 0. _d 0            cg2d_r(i,j,bi,bj) = 0. _d 0
73            cg2d_s(I,J,bi,bj) = 0. _d 0            cg2d_s(i,j,bi,bj) = 0. _d 0
74  #ifdef ALLOW_CG2D_NSA  #ifdef ALLOW_CG2D_NSA
75            cg2d_z(I,J,bi,bj) = 0. _d 0            cg2d_z(i,j,bi,bj) = 0. _d 0
76  #endif /* ALLOW_CG2D_NSA */  #endif /* ALLOW_CG2D_NSA */
77  #ifdef ALLOW_SRCG  #ifdef ALLOW_SRCG
78            cg2d_y(I,J,bi,bj) = 0. _d 0            cg2d_y(i,j,bi,bj) = 0. _d 0
79            cg2d_v(I,J,bi,bj) = 0. _d 0            cg2d_v(i,j,bi,bj) = 0. _d 0
80  #endif /*  ALLOW_SRCG */  #endif /*  ALLOW_SRCG */
81           ENDDO           ENDDO
82          ENDDO          ENDDO
# Line 89  C     aS2d: integral in Z Ay/dY Line 89  C     aS2d: integral in Z Ay/dY
89        myNorm = 0. _d 0        myNorm = 0. _d 0
90        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
91         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
92          DO J=1,sNy          DO j=1,sNy
93           DO I=1,sNx           DO i=1,sNx
94            aW2d(I,J,bi,bj) = 0. _d 0            aW2d(i,j,bi,bj) = 0. _d 0
95            aS2d(I,J,bi,bj) = 0. _d 0            aS2d(i,j,bi,bj) = 0. _d 0
96           ENDDO           ENDDO
97          ENDDO          ENDDO
98          DO K=1,Nr          DO k=1,Nr
99           DO J=1,sNy           DO j=1,sNy
100            DO I=1,sNx            DO i=1,sNx
101  C  deep-model: *deepFacC (faceArea), /deepFacC (recip_dx,y): => no net effect  C  deep-model: *deepFacC (faceArea), /deepFacC (recip_dx,y): => no net effect
102             faceArea = _dyG(I,J,bi,bj)*drF(K)             faceArea = _dyG(i,j,bi,bj)*drF(k)
103       &               *_hFacW(I,J,K,bi,bj)       &               *_hFacW(i,j,k,bi,bj)
104             aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj)             aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)
105       &              + implicSurfPress*implicDiv2DFlow       &              + implicSurfPress*implicDiv2DFlow
106       &               *faceArea*recip_dxC(I,J,bi,bj)       &               *faceArea*recip_dxC(i,j,bi,bj)
107             faceArea = _dxG(I,J,bi,bj)*drF(K)             faceArea = _dxG(i,j,bi,bj)*drF(k)
108       &               *_hFacS(I,J,K,bi,bj)       &               *_hFacS(i,j,k,bi,bj)
109             aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj)             aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)
110       &              + implicSurfPress*implicDiv2DFlow       &              + implicSurfPress*implicDiv2DFlow
111       &               *faceArea*recip_dyC(I,J,bi,bj)       &               *faceArea*recip_dyC(i,j,bi,bj)
112            ENDDO            ENDDO
113           ENDDO           ENDDO
114          ENDDO          ENDDO
115  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
116          IF (useOBCS) THEN          IF (useOBCS) THEN
117           DO I=1,sNx           DO i=1,sNx
118            IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj),bi,bj)=0.            IF (OB_Jn(i,bi,bj).NE.0) aS2d(i, OB_Jn(i,bi,bj), bi,bj)=0.
119            IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj)+1,bi,bj)=0.            IF (OB_Jn(i,bi,bj).NE.0) aS2d(i,OB_Jn(i,bi,bj)+1,bi,bj)=0.
120            IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj)+1,bi,bj)=0.            IF (OB_Js(i,bi,bj).NE.0) aS2d(i,OB_Js(i,bi,bj)+1,bi,bj)=0.
121            IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj),bi,bj)=0.            IF (OB_Js(i,bi,bj).NE.0) aS2d(i, OB_Js(i,bi,bj), bi,bj)=0.
122           ENDDO           ENDDO
123           DO J=1,sNy           DO j=1,sNy
124            IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj),J,bi,bj)=0.            IF (OB_Ie(j,bi,bj).NE.0) aW2d(OB_Ie(j,bi,bj), j, bi,bj)=0.
125            IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj)+1,J,bi,bj)=0.            IF (OB_Ie(j,bi,bj).NE.0) aW2d(OB_Ie(j,bi,bj)+1,j,bi,bj)=0.
126            IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj)+1,J,bi,bj)=0.            IF (OB_Iw(j,bi,bj).NE.0) aW2d(OB_Iw(j,bi,bj)+1,j,bi,bj)=0.
127            IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj),J,bi,bj)=0.            IF (OB_Iw(j,bi,bj).NE.0) aW2d(OB_Iw(j,bi,bj), j, bi,bj)=0.
128           ENDDO           ENDDO
129          ENDIF          ENDIF
130  #endif  #endif
131          DO J=1,sNy          DO j=1,sNy
132           DO I=1,sNx           DO i=1,sNx
133            myNorm = MAX(ABS(aW2d(I,J,bi,bj)),myNorm)            myNorm = MAX(ABS(aW2d(i,j,bi,bj)),myNorm)
134            myNorm = MAX(ABS(aS2d(I,J,bi,bj)),myNorm)            myNorm = MAX(ABS(aS2d(i,j,bi,bj)),myNorm)
135           ENDDO           ENDDO
136          ENDDO          ENDDO
137         ENDDO         ENDDO
# Line 144  C  deep-model: *deepFacC (faceArea), /de Line 144  C  deep-model: *deepFacC (faceArea), /de
144        ENDIF        ENDIF
145        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
146         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
147          DO J=1,sNy          DO j=1,sNy
148           DO I=1,sNx           DO i=1,sNx
149            aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj)*myNorm            aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)*myNorm
150            aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj)*myNorm            aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)*myNorm
151           ENDDO           ENDDO
152          ENDDO          ENDDO
153         ENDDO         ENDDO
# Line 210  C           For now I have introduced a Line 210  C           For now I have introduced a
210  C           defaults to 0.51 but can be set at runtime.  C           defaults to 0.51 but can be set at runtime.
211        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
212         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
213          DO J=1,sNy  C-    calculate and store solver main diagonal :
214           DO I=1,sNx          DO j=0,sNy+1
215            ks = ksurfC(I,J,bi,bj)           DO i=0,sNx+1
216            pC(I,J,bi,bj) = 1. _d 0             ks = ksurfC(i,j,bi,bj)
217            aC = -(             aC2d(i,j,bi,bj) = -(
218       &     aW2d(I,J,bi,bj) + aW2d(I+1,J  ,bi,bj)       &       aW2d(i,j,bi,bj) + aW2d(i+1,j, bi,bj)
219       &    +aS2d(I,J,bi,bj) + aS2d(I  ,J+1,bi,bj)       &      +aS2d(i,j,bi,bj) + aS2d( i,j+1,bi,bj)
220       &    +freeSurfFac*myNorm*recip_Bo(I,J,bi,bj)*deepFac2F(ks)       &      +freeSurfFac*myNorm*recip_Bo(i,j,bi,bj)*deepFac2F(ks)
221       &                *rA(I,J,bi,bj)/deltaTMom/deltaTfreesurf       &                  *rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
222       &    )       &                        )
223            aCs = -(           ENDDO
224       &     aW2d(I,J-1,bi,bj) + aW2d(I+1,J-1,bi,bj)          ENDDO
225       &    +aS2d(I,J-1,bi,bj) + aS2d(I  ,J  ,bi,bj)          DO j=1,sNy
226       &    +freeSurfFac*myNorm*recip_Bo(I,J-1,bi,bj)*deepFac2F(ks)           DO i=1,sNx
227       &                *rA(I,J-1,bi,bj)/deltaTMom/deltaTfreesurf             aC  = aC2d( i, j, bi,bj)
228       &    )             aCs = aC2d( i,j-1,bi,bj)
229            aCw = -(             aCw = aC2d(i-1,j, bi,bj)
230       &     aW2d(I-1,J,bi,bj) + aW2d(I  ,J  ,bi,bj)             IF ( aC .EQ. 0. ) THEN
231       &    +aS2d(I-1,J,bi,bj) + aS2d(I-1,J+1,bi,bj)              pC(i,j,bi,bj) = 1. _d 0
232       &    +freeSurfFac*myNorm*recip_Bo(I-1,J,bi,bj)*deepFac2F(ks)             ELSE
233       &                *rA(I-1,J,bi,bj)/deltaTMom/deltaTfreesurf              pC(i,j,bi,bj) =  1. _d 0 / aC
234       &    )             ENDIF
235            IF ( aC .EQ. 0. ) THEN             IF ( aC + aCw .EQ. 0. ) THEN
236              pC(I,J,bi,bj) = 1. _d 0              pW(i,j,bi,bj) = 0.
237            ELSE             ELSE
238             pC(I,J,bi,bj) =  1. _d 0 / aC              pW(i,j,bi,bj) =
239            ENDIF       &         -aW2d(i,j,bi,bj)/((cg2dpcOffDFac *(aCw+aC))**2 )
240            IF ( aC + aCw .EQ. 0. ) THEN             ENDIF
241             pW(I,J,bi,bj) = 0.             IF ( aC + aCs .EQ. 0. ) THEN
242            ELSE              pS(i,j,bi,bj) = 0.
243             pW(I,J,bi,bj) =             ELSE
244       &     -aW2d(I  ,J  ,bi,bj)/((cg2dpcOffDFac *(aCw+aC))**2 )              pS(i,j,bi,bj) =
245            ENDIF       &         -aS2d(i,j,bi,bj)/((cg2dpcOffDFac *(aCs+aC))**2 )
246            IF ( aC + aCs .EQ. 0. ) THEN             ENDIF
247             pS(I,J,bi,bj) = 0.  C          pC(i,j,bi,bj) = 1.
248            ELSE  C          pW(i,j,bi,bj) = 0.
249             pS(I,J,bi,bj) =  C          pS(i,j,bi,bj) = 0.
      &     -aS2d(I  ,J  ,bi,bj)/((cg2dpcOffDFac *(aCs+aC))**2 )  
           ENDIF  
 C-    store solver main diagonal :  
           aC2d(I,J,bi,bj) = aC  
 C         pC(I,J,bi,bj) = 1.  
 C         pW(I,J,bi,bj) = 0.  
 C         pS(I,J,bi,bj) = 0.  
250           ENDDO           ENDDO
251          ENDDO          ENDDO
252         ENDDO         ENDDO
Legend:
Removed from v.1.48  
changed lines
  Added in v.1.49
  ViewVC Help
Powered by ViewVC 1.1.22