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

Diff of /MITgcm/model/src/ini_spherical_polar_grid.F

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

revision 1.6 by adcroft, Mon Jun 22 15:26:25 1998 UTC revision 1.11 by adcroft, Mon Nov 30 23:45:25 1998 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    
3  #include "CPP_EEOPTIONS.h"  #include "CPP_OPTIONS.h"
4    
5  CStartOfInterface  CStartOfInterface
6        SUBROUTINE INI_SPHERICAL_POLAR_GRID( myThid )        SUBROUTINE INI_SPHERICAL_POLAR_GRID( myThid )
# Line 65  C     latS         values. Line 65  C     latS         values.
65  C     I,J,K  C     I,J,K
66        _RL    xG, yG, zG        _RL    xG, yG, zG
67        _RL    phi        _RL    phi
68        _RL    zUpper(Nz), zLower(Nz)        _RL    zUpper(Nr), zLower(Nr)
69        _RL    xBase, yBase        _RL    xBase, yBase
70        INTEGER iG, jG        INTEGER iG, jG
71        INTEGER bi, bj        INTEGER bi, bj
# Line 102  C           longitude - phiMin and theta Line 102  C           longitude - phiMin and theta
102            xc(I,J,bi,bj)  = xG + delX(iG+i-1)*0.5 _d 0            xc(I,J,bi,bj)  = xG + delX(iG+i-1)*0.5 _d 0
103            yc(I,J,bi,bj)  = yG + delY(jG+j-1)*0.5 _d 0            yc(I,J,bi,bj)  = yG + delY(jG+j-1)*0.5 _d 0
104            xG = xG + delX(iG+I-1)            xG = xG + delX(iG+I-1)
105            dxF(I,J,bi,bj) = delX(iG+i-1)*deg2rad*rSphere*COS(yc(I,J,bi,bj)*deg2rad)            dxF(I,J,bi,bj) = delX(iG+i-1)*deg2rad
106         &    *rSphere*COS(yc(I,J,bi,bj)*deg2rad)
107            dyF(I,J,bi,bj) = delY(jG+j-1)*deg2rad*rSphere            dyF(I,J,bi,bj) = delY(jG+j-1)*deg2rad*rSphere
108           ENDDO           ENDDO
109           yG = yG + delY(jG+J-1)           yG = yG + delY(jG+J-1)
# Line 134  C     dxG, dyG are separations between c Line 135  C     dxG, dyG are separations between c
135        ENDDO        ENDDO
136        _EXCH_XY_R4(dxG, myThid )        _EXCH_XY_R4(dxG, myThid )
137        _EXCH_XY_R4(dyG, myThid )        _EXCH_XY_R4(dyG, myThid )
 C     dxV, dyU are separations between velocity points along cell faces.  
       DO bj = myByLo(myThid), myByHi(myThid)  
        DO bi = myBxLo(myThid), myBxHi(myThid)  
         DO J=1,sNy  
          DO I=1,sNx  
           dxV(I,J,bi,bj) = (dxG(I,J,bi,bj)+dxG(I-1,J,bi,bj))*0.5 _d 0  
           dyU(I,J,bi,bj) = (dyG(I,J,bi,bj)+dyG(I,J-1,bi,bj))*0.5 _d 0  
          ENDDO  
         ENDDO  
        ENDDO  
       ENDDO  
       _EXCH_XY_R4(dxV, myThid )  
       _EXCH_XY_R4(dyU, myThid )  
138  C     dxC, dyC is separation between cell centers  C     dxC, dyC is separation between cell centers
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)
# Line 160  C     dxC, dyC is separation between cel Line 148  C     dxC, dyC is separation between cel
148        ENDDO        ENDDO
149        _EXCH_XY_R4(dxC, myThid )        _EXCH_XY_R4(dxC, myThid )
150        _EXCH_XY_R4(dyC, myThid )        _EXCH_XY_R4(dyC, myThid )
151  C     Calculate recipricols  C     dxV, dyU are separations between velocity points along cell faces.
152        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
153         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
154          DO J=1,sNy          DO J=1,sNy
155           DO I=1,sNx           DO I=1,sNx
156            rDxG(I,J,bi,bj)=1.d0/dxG(I,J,bi,bj)            dxV(I,J,bi,bj) = (dxG(I,J,bi,bj)+dxG(I-1,J,bi,bj))*0.5 _d 0
157            rDyG(I,J,bi,bj)=1.d0/dyG(I,J,bi,bj)  #ifdef OLD_UV_GEOMETRY
158            rDxC(I,J,bi,bj)=1.d0/dxC(I,J,bi,bj)            dyU(I,J,bi,bj) = (dyG(I,J,bi,bj)+dyG(I,J-1,bi,bj))*0.5 _d 0
159            rDyC(I,J,bi,bj)=1.d0/dyC(I,J,bi,bj)  #else
160            rDxF(I,J,bi,bj)=1.d0/dxF(I,J,bi,bj)            dyU(I,J,bi,bj) = (dyC(I,J,bi,bj)+dyC(I-1,J,bi,bj))*0.5 _d 0
161            rDyF(I,J,bi,bj)=1.d0/dyF(I,J,bi,bj)  #endif
           rDxV(I,J,bi,bj)=1.d0/dxV(I,J,bi,bj)  
           rDyU(I,J,bi,bj)=1.d0/dyU(I,J,bi,bj)  
162           ENDDO           ENDDO
163          ENDDO          ENDDO
164         ENDDO         ENDDO
165        ENDDO        ENDDO
166        _EXCH_XY_R4(rDxG, myThid )        _EXCH_XY_R4(dxV, myThid )
167        _EXCH_XY_R4(rDyG, myThid )        _EXCH_XY_R4(dyU, myThid )
       _EXCH_XY_R4(rDxC, myThid )  
       _EXCH_XY_R4(rDyC, myThid )  
       _EXCH_XY_R4(rDxF, myThid )  
       _EXCH_XY_R4(rDyF, myThid )  
       _EXCH_XY_R4(rDxV, myThid )  
       _EXCH_XY_R4(rDyU, myThid )  
168  C     Calculate vertical face area and trigonometric terms  C     Calculate vertical face area and trigonometric terms
169        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
170         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
# Line 193  C     Calculate vertical face area and t Line 173  C     Calculate vertical face area and t
173            jG = myYGlobalLo + (bj-1)*sNy + J-1            jG = myYGlobalLo + (bj-1)*sNy + J-1
174            latS = yc(i,j,bi,bj)-delY(jG)*0.5 _d 0            latS = yc(i,j,bi,bj)-delY(jG)*0.5 _d 0
175            latN = yc(i,j,bi,bj)+delY(jG)*0.5 _d 0            latN = yc(i,j,bi,bj)+delY(jG)*0.5 _d 0
176            zA(I,J,bi,bj) = dyF(I,J,bi,bj)  #ifdef OLD_UV_GEOMETRY
177              rA(I,J,bi,bj) = dyF(I,J,bi,bj)
178         &    *rSphere*(SIN(latN*deg2rad)-SIN(latS*deg2rad))
179    #else
180              rA(I,J,bi,bj) = rSphere*delX(iG)*deg2rad
181       &    *rSphere*(SIN(latN*deg2rad)-SIN(latS*deg2rad))       &    *rSphere*(SIN(latN*deg2rad)-SIN(latS*deg2rad))
182    #endif
183    C         Area cannot be zero but sin can be if lat if < -90.
184              IF ( rA(I,J,bi,bj) .LT. 0. ) rA(I,J,bi,bj) = -rA(I,J,bi,bj)
185            tanPhiAtU(i,j,bi,bj)=tan(_yC(i,j,bi,bj)*deg2rad)            tanPhiAtU(i,j,bi,bj)=tan(_yC(i,j,bi,bj)*deg2rad)
186            tanPhiAtV(i,j,bi,bj)=tan(latS*deg2rad)            tanPhiAtV(i,j,bi,bj)=tan(latS*deg2rad)
187           ENDDO           ENDDO
188          ENDDO          ENDDO
189         ENDDO         ENDDO
190        ENDDO        ENDDO
191          _EXCH_XY_R4 (rA       , myThid )
192          _EXCH_XY_R4 (tanPhiAtU , myThid )
193          _EXCH_XY_R4 (tanPhiAtV , myThid )
194        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
195         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
196          DO K=1,Nz          DO J=1,sNy
197           DO J=1,sNy           DO I=1,sNx
198            DO I=1,sNx            iG = myXGlobalLo + (bi-1)*sNx + I-1
199             IF (HFacC(I,J,K,bi,bj) .NE. 0. D0 ) THEN            jG = myYGlobalLo + (bj-1)*sNy + J-1
200              rHFacC(I,J,K,bi,bj) = 1. D0 / HFacC(I,J,K,bi,bj)            latS = yc(i,j-1,bi,bj)
201             ELSE            latN = yc(i,j,bi,bj)
202              rHFacC(I,J,K,bi,bj) = 0. D0  #ifdef OLD_UV_GEOMETRY
203             ENDIF            rAw(I,J,bi,bj) = 0.5*(rA(I,J,bi,bj)+rA(I-1,J,bi,bj))
204             IF (HFacW(I,J,K,bi,bj) .NE. 0. D0 ) THEN            rAs(I,J,bi,bj) = 0.5*(rA(I,J,bi,bj)+rA(I,J-1,bi,bj))
205              rHFacW(I,J,K,bi,bj) = 1. D0 / HFacW(I,J,K,bi,bj)  #else
206              maskW(I,J,K,bi,bj) = 1. D0            rAw(I,J,bi,bj) = 0.5*(rA(I,J,bi,bj)+rA(I-1,J,bi,bj))
207             ELSE            rAs(I,J,bi,bj) = rSphere*delX(iG)*deg2rad
208              rHFacW(I,J,K,bi,bj) = 0. D0       &    *rSphere*(SIN(latN*deg2rad)-SIN(latS*deg2rad))
209              maskW(I,J,K,bi,bj) = 0.0 D0  #endif
            ENDIF  
            IF (HFacS(I,J,K,bi,bj) .NE. 0. D0 ) THEN  
             rHFacS(I,J,K,bi,bj) = 1. D0 / HFacS(I,J,K,bi,bj)  
             maskS(I,J,K,bi,bj) = 1. D0  
            ELSE  
             rHFacS(I,J,K,bi,bj) = 0. D0  
             maskS(I,J,K,bi,bj) = 0. D0  
            ENDIF  
           ENDDO  
210           ENDDO           ENDDO
211          ENDDO          ENDDO
212         ENDDO         ENDDO
213        ENDDO        ENDDO
214  C     Now sync. and get/send edge regions that are shared with        _EXCH_XY_R4 (rAw      , myThid )
215  C     other threads.        _EXCH_XY_R4 (rAs      , myThid )
       _EXCH_XYZ_R4(rHFacC    , myThid )  
       _EXCH_XYZ_R4(rHFacW    , myThid )  
       _EXCH_XYZ_R4(rHFacS    , myThid )  
       _EXCH_XYZ_R4(maskW    , myThid )  
       _EXCH_XYZ_R4(maskS    , myThid )  
       _EXCH_XY_R4 (zA       , myThid )  
       _EXCH_XY_R4 (tanPhiAtU , myThid )  
       _EXCH_XY_R4 (tanPhiAtV , myThid )  
   
216  C  C
217        RETURN        RETURN
218        END        END
Legend:
Removed from v.1.6  
changed lines
  Added in v.1.11
  ViewVC Help
Powered by ViewVC 1.1.22