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

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

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

revision 1.17 by adcroft, Tue May 18 17:57:01 1999 UTC revision 1.44 by jmc, Thu Dec 17 23:47:06 2009 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4    #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    
7  CStartOfInterface  CBOP
8    C     !ROUTINE: INI_MASKS_ETC
9    C     !INTERFACE:
10        SUBROUTINE INI_MASKS_ETC( myThid )        SUBROUTINE INI_MASKS_ETC( myThid )
11  C     /==========================================================\  C     !DESCRIPTION: \bv
12  C     | SUBROUTINE INI_MASKS_ETC                                 |  C     *==========================================================*
13  C     | o Initialise masks and topography factors                |  C     | SUBROUTINE INI_MASKS_ETC
14  C     |==========================================================|  C     | o Initialise masks and topography factors
15  C     | These arrays are used throughout the code and describe   |  C     *==========================================================*
16  C     | the topography of the domain through masks (0s and 1s)   |  C     | These arrays are used throughout the code and describe
17  C     | and fractional height factors (0<hFac<1). The latter     |  C     | the topography of the domain through masks (0s and 1s)
18  C     | distinguish between the lopped-cell and full-step        |  C     | and fractional height factors (0<hFac<1). The latter
19  C     | topographic representations.                             |  C     | distinguish between the lopped-cell and full-step
20  C     \==========================================================/  C     | topographic representations.
21        IMPLICIT NONE  C     *==========================================================*
22    C     \ev
23    
24    C     !USES:
25          IMPLICIT NONE
26  C     === Global variables ===  C     === Global variables ===
27  #include "SIZE.h"  #include "SIZE.h"
28  #include "EEPARAMS.h"  #include "EEPARAMS.h"
29  #include "PARAMS.h"  #include "PARAMS.h"
30  #include "GRID.h"  #include "GRID.h"
31    #include "SURFACE.h"
32    
33    C     !INPUT/OUTPUT PARAMETERS:
34  C     == Routine arguments ==  C     == Routine arguments ==
35  C     myThid -  Number of this instance of INI_MASKS_ETC  C     myThid  ::  Number of this instance of INI_MASKS_ETC
36        INTEGER myThid        INTEGER myThid
 CEndOfInterface  
37    
38    C     !LOCAL VARIABLES:
39  C     == Local variables ==  C     == Local variables ==
40  C     bi,bj  - Loop counters  C     bi,bj   :: tile indices
41  C     I,J,K  C     I,J,K   :: Loop counters
42    C     tmpfld  :: Temporary array used to compute & write Total Depth
43          _RS tmpfld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
44        INTEGER bi, bj        INTEGER bi, bj
45        INTEGER  I, J, K        INTEGER I, J, K
46  #ifdef ALLOW_NONHYDROSTATIC        _RL hFacCtmp
47        INTEGER Km1        _RL hFacMnSz
48        _RL hFacUpper,hFacLower  CEOP
 #endif  
49    
50  C     Calculate lopping factor hFacC  C- Calculate lopping factor hFacC : over-estimate the part inside of the domain
51    C    taking into account the lower_R Boundary (Bathymetrie / Top of Atmos)
52        DO bj=myByLo(myThid), myByHi(myThid)        DO bj=myByLo(myThid), myByHi(myThid)
53         DO bi=myBxLo(myThid), myBxHi(myThid)         DO bi=myBxLo(myThid), myBxHi(myThid)
54          DO K=1, Nr          DO K=1, Nr
55           DO J=1,sNy           hFacMnSz=max( hFacMin, min(hFacMinDr*recip_drF(k),1. _d 0) )
56            DO I=1,sNx           DO J=1-Oly,sNy+Oly
57  C          Round depths within a small fraction of layer depth to that            DO I=1-Olx,sNx+Olx
58  C          layer depth.  C      o Non-dimensional distance between grid bound. and domain lower_R bound.
59             IF ( ABS(H(I,J,bi,bj)-rF(K)) .LT.             hFacCtmp = (rF(K)-R_low(I,J,bi,bj))*recip_drF(K)
60       &          1. _d -6*ABS(rF(K)) .AND.  C      o Select between, closed, open or partial (0,1,0-1)
61       &          ABS(H(I,J,bi,bj)-rF(K)) .LT.             hFacCtmp=min( max( hFacCtmp, 0. _d 0) , 1. _d 0)
62       &          1. _d -6*ABS(H(I,J,bi,bj)) )THEN  C      o Impose minimum fraction and/or size (dimensional)
63              H(I,J,bi,bj) = rF(K)             IF (hFacCtmp.LT.hFacMnSz) THEN
64             ENDIF              IF (hFacCtmp.LT.hFacMnSz*0.5) THEN
            IF     ( H(I,J,bi,bj)*rkFac .GE. rF(K)*rkFac ) THEN  
 C           Top of cell is below base of domain  
             hFacC(I,J,K,bi,bj) = 0.  
            ELSEIF ( H(I,J,bi,bj)*rkFac .LE. rF(K+1)*rkFac ) THEN  
 C           Base of domain is below bottom of this cell  
             hFacC(I,J,K,bi,bj) = 1.  
            ELSE  
 C           Base of domain is in this cell  
 C           Set hFac to the fraction of the cell that is open.  
             hFacC(I,J,K,bi,bj) =  
      &       (rF(K)*rkFac-H(I,J,bi,bj)*rkFac)*recip_drF(K)  
            ENDIF  
 C Impose minimum fraction  
            IF (hFacC(I,J,K,bi,bj).LT.hFacMin) THEN  
             IF (hFacC(I,J,K,bi,bj).LT.hFacMin*0.5) THEN  
65               hFacC(I,J,K,bi,bj)=0.               hFacC(I,J,K,bi,bj)=0.
66              ELSE              ELSE
67               hFacC(I,J,K,bi,bj)=hFacMin               hFacC(I,J,K,bi,bj)=hFacMnSz
68              ENDIF              ENDIF
69               ELSE
70                 hFacC(I,J,K,bi,bj)=hFacCtmp
71             ENDIF             ENDIF
72  C Impose minimum size (dimensional)            ENDDO
73             IF (drF(k)*hFacC(I,J,K,bi,bj).LT.hFacMinDr) THEN           ENDDO
74              IF (drF(k)*hFacC(I,J,K,bi,bj).LT.hFacMinDr*0.5) THEN          ENDDO
75    
76    C-  Re-calculate lower-R Boundary position, taking into account hFacC
77            DO J=1-Oly,sNy+Oly
78             DO I=1-Olx,sNx+Olx
79              R_low(I,J,bi,bj) = rF(1)
80              DO K=Nr,1,-1
81               R_low(I,J,bi,bj) = R_low(I,J,bi,bj)
82         &                      - drF(k)*hFacC(I,J,K,bi,bj)
83              ENDDO
84             ENDDO
85            ENDDO
86    C - end bi,bj loops.
87           ENDDO
88          ENDDO
89    
90    C-  Calculate lopping factor hFacC : Remove part outside of the domain
91    C    taking into account the Reference (=at rest) Surface Position Ro_surf
92          DO bj=myByLo(myThid), myByHi(myThid)
93           DO bi=myBxLo(myThid), myBxHi(myThid)
94            DO K=1, Nr
95             hFacMnSz=max( hFacMin, min(hFacMinDr*recip_drF(k),1. _d 0) )
96             DO J=1-Oly,sNy+Oly
97              DO I=1-Olx,sNx+Olx
98    C      o Non-dimensional distance between grid boundary and model surface
99               hFacCtmp = (rF(k)-Ro_surf(I,J,bi,bj))*recip_drF(K)
100    C      o Reduce the previous fraction : substract the outside part.
101               hFacCtmp = hFacC(I,J,K,bi,bj) - max( hFacCtmp, 0. _d 0)
102    C      o set to zero if empty Column :
103               hFacCtmp = max( hFacCtmp, 0. _d 0)
104    C      o Impose minimum fraction and/or size (dimensional)
105               IF (hFacCtmp.LT.hFacMnSz) THEN
106                IF (hFacCtmp.LT.hFacMnSz*0.5) THEN
107               hFacC(I,J,K,bi,bj)=0.               hFacC(I,J,K,bi,bj)=0.
108              ELSE              ELSE
109               hFacC(I,J,K,bi,bj)=hFacMinDr*recip_drF(k)               hFacC(I,J,K,bi,bj)=hFacMnSz
110              ENDIF              ENDIF
111               ELSE
112                 hFacC(I,J,K,bi,bj)=hFacCtmp
113             ENDIF             ENDIF
            depthInK(i,j,bi,bj) = depthInK(i,j,bi,bj) + 1.  
 Crg  &                          +hFacC(i,j,k,bi,bj)  
114            ENDDO            ENDDO
115           ENDDO           ENDDO
116          ENDDO          ENDDO
117         ENDDO         ENDDO
118        ENDDO        ENDDO
       _EXCH_XYZ_R4(hFacC , myThid )  
       _EXCH_XY_R4( depthInK, myThid )  
119    
120        CALL PLOT_FIELD_XYRS( depthInK,  #ifdef ALLOW_SHELFICE
121       & 'Model Depths K Index' , 1, myThid )        IF ( useShelfIce ) THEN
122    C--   Modify lopping factor hFacC : Remove part outside of the domain
123  C Re-calculate depth of ocean, taking into account hFacC  C     taking into account the Reference (=at rest) Surface Position Ro_shelfIce
124           CALL SHELFICE_UPDATE_MASKS(
125         I     rF, recip_drF,
126         U     hFacC,
127         I     myThid )
128          ENDIF
129    #endif /* ALLOW_SHELFICE */
130    
131    C-  Re-calculate Reference surface position, taking into account hFacC
132    C   initialize Total column fluid thickness and surface k index
133    C       Note: if no fluid (continent) ==> ksurf = Nr+1
134        DO bj=myByLo(myThid), myByHi(myThid)        DO bj=myByLo(myThid), myByHi(myThid)
135         DO bi=myBxLo(myThid), myBxHi(myThid)         DO bi=myBxLo(myThid), myBxHi(myThid)
136          DO J=1,sNy          DO J=1-Oly,sNy+Oly
137           DO I=1,sNx           DO I=1-Olx,sNx+Olx
138            H(I,J,bi,bj)=0.            tmpfld(I,J,bi,bj) = 0.
139            DO K=1,Nr            ksurfC(I,J,bi,bj) = Nr+1
140             H(I,J,bi,bj)=H(I,J,bi,bj)-            maskH(I,J,bi,bj)  = 0.
141       &       rkFac*drF(k)*hFacC(I,J,K,bi,bj)            Ro_surf(I,J,bi,bj) = R_low(I,J,bi,bj)
142              DO K=Nr,1,-1
143               Ro_surf(I,J,bi,bj) = Ro_surf(I,J,bi,bj)
144         &                        + drF(k)*hFacC(I,J,K,bi,bj)
145               IF (hFacC(I,J,K,bi,bj).NE.0.) THEN
146                ksurfC(I,J,bi,bj) = k
147                maskH(I,J,bi,bj)  = 1.
148                tmpfld(I,J,bi,bj) = tmpfld(I,J,bi,bj) + 1.
149               ENDIF
150              ENDDO
151              kLowC(I,J,bi,bj) = 0
152              DO K= 1, Nr
153               IF (hFacC(I,J,K,bi,bj).NE.0) THEN
154                  kLowC(I,J,bi,bj) = K
155               ENDIF
156            ENDDO            ENDDO
157              maskInC(I,J,bi,bj)= 0.
158              IF ( ksurfC(I,J,bi,bj).LE.Nr ) maskInC(I,J,bi,bj)= 1.
159           ENDDO           ENDDO
160          ENDDO          ENDDO
161    C - end bi,bj loops.
162         ENDDO         ENDDO
163        ENDDO        ENDDO
164        _EXCH_XY_R4(H , myThid )  
165        CALL WRITE_FLD_XY_RS( 'Depth',' ',H,0,myThid)  C     CALL PLOT_FIELD_XYRS( tmpfld,
166  C     CALL MDSWRITEFIELD( 'Depth', writeBinaryPrec, .TRUE.,  C    &         'Model Depths K Index' , 1, myThid )
167  C    &                    'RS', 1, H, 1, -1, myThid )        CALL PLOT_FIELD_XYRS(R_low,
168         &         'Model R_low (ini_masks_etc)', 1, myThid)
169          CALL PLOT_FIELD_XYRS(Ro_surf,
170         &         'Model Ro_surf (ini_masks_etc)', 1, myThid)
171    
172  C     Calculate quantities derived from XY depth map  C     Calculate quantities derived from XY depth map
173        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
174         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
175          DO J=1,sNy          DO j=1-Oly,sNy+Oly
176           DO I=1,sNx           DO i=1-Olx,sNx+Olx
177  C         Inverse of depth  C         Total fluid column thickness (r_unit) :
178            IF ( h(i,j,bi,bj) .EQ. 0. _d 0 ) THEN  c         Rcolumn(i,j,bi,bj)= Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
179             recip_H(i,j,bi,bj) = 0. _d 0            tmpfld(i,j,bi,bj) = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
180    C         Inverse of fluid column thickness (1/r_unit)
181              IF ( tmpfld(i,j,bi,bj) .LE. 0. ) THEN
182               recip_Rcol(i,j,bi,bj) = 0.
183            ELSE            ELSE
184             recip_H(i,j,bi,bj) = 1. _d 0 /  abs( H(i,j,bi,bj) )             recip_Rcol(i,j,bi,bj) = 1. _d 0 / tmpfld(i,j,bi,bj)
185            ENDIF            ENDIF
           depthInK(i,j,bi,bj) = 0.  
186           ENDDO           ENDDO
187          ENDDO          ENDDO
188         ENDDO         ENDDO
189        ENDDO        ENDDO
190        _EXCH_XY_R4(   recip_H, myThid )  C     _EXCH_XY_RS(   recip_Rcol, myThid )
191    
192  C     hFacW and hFacS (at U and V points)  C     hFacW and hFacS (at U and V points)
193        DO bj=myByLo(myThid), myByHi(myThid)        DO bj=myByLo(myThid), myByHi(myThid)
194         DO bi=myBxLo(myThid), myBxHi(myThid)         DO bi=myBxLo(myThid), myBxHi(myThid)
195          DO K=1, Nr          DO K=1, Nr
196           DO J=1,sNy           DO J=1-Oly,sNy+Oly
197            DO I=1,sNx            hFacW(1-OLx,J,K,bi,bj)= 0.
198              DO I=2-Olx,sNx+Olx
199             hFacW(I,J,K,bi,bj)=             hFacW(I,J,K,bi,bj)=
200       &       MIN(hFacC(I,J,K,bi,bj),hFacC(I-1,J,K,bi,bj))       &       MIN(hFacC(I,J,K,bi,bj),hFacC(I-1,J,K,bi,bj))
201              ENDDO
202             ENDDO
203             DO I=1-Olx,sNx+Olx
204               hFacS(I,1-OLy,K,bi,bj)= 0.
205             ENDDO
206             DO J=2-Oly,sNy+oly
207              DO I=1-Olx,sNx+Olx
208             hFacS(I,J,K,bi,bj)=             hFacS(I,J,K,bi,bj)=
209       &       MIN(hFacC(I,J,K,bi,bj),hFacC(I,J-1,K,bi,bj))       &       MIN(hFacC(I,J,K,bi,bj),hFacC(I,J-1,K,bi,bj))
210            ENDDO            ENDDO
# Line 144  C     hFacW and hFacS (at U and V points Line 212  C     hFacW and hFacS (at U and V points
212          ENDDO          ENDDO
213         ENDDO         ENDDO
214        ENDDO        ENDDO
215        _EXCH_XYZ_R4(hFacW , myThid )        CALL EXCH_UV_XYZ_RS(hFacW,hFacS,.FALSE.,myThid)
216        _EXCH_XYZ_R4(hFacS , myThid )  C     The following block allows thin walls representation of non-periodic
217    C     boundaries such as happen on the lat-lon grid at the N/S poles.
218    C     We should really supply a flag for doing this.
219          DO bj=myByLo(myThid), myByHi(myThid)
220           DO bi=myBxLo(myThid), myBxHi(myThid)
221            DO K=1, Nr
222             DO J=1-Oly,sNy+Oly
223              DO I=1-Olx,sNx+Olx
224               IF (dyG(I,J,bi,bj).EQ.0.) hFacW(I,J,K,bi,bj)=0.
225               IF (dxG(I,J,bi,bj).EQ.0.) hFacS(I,J,K,bi,bj)=0.
226              ENDDO
227             ENDDO
228            ENDDO
229           ENDDO
230          ENDDO
231    
232    #ifdef ALLOW_NONHYDROSTATIC
233    C     rLow & reference rSurf at Western & Southern edges (U and V points)
234          DO bj=myByLo(myThid), myByHi(myThid)
235           DO bi=myBxLo(myThid), myBxHi(myThid)
236            I = 1-OlX
237            DO J=1-Oly,sNy+Oly
238               rLowW (i,j,bi,bj) = 0.
239               rSurfW(i,j,bi,bj) = 0.
240            ENDDO
241            J = 1-Oly
242            DO I=1-Olx,sNx+Olx
243               rLowS (i,j,bi,bj) = 0.
244               rSurfS(i,j,bi,bj) = 0.
245            ENDDO
246            DO J=1-Oly,sNy+Oly
247             DO I=2-Olx,sNx+Olx
248               rSurfW(i,j,bi,bj) =
249         &           MIN( Ro_surf(i-1,j,bi,bj), Ro_surf(i,j,bi,bj) )
250               rLowW(i,j,bi,bj)  =
251         &           MAX(   R_low(i-1,j,bi,bj),   R_low(i,j,bi,bj) )
252             ENDDO
253            ENDDO
254            DO J=2-Oly,sNy+Oly
255             DO I=1-Olx,sNx+Olx
256               rSurfS(i,j,bi,bj) =
257         &           MIN( Ro_surf(i,j-1,bi,bj), Ro_surf(i,j,bi,bj) )
258               rLowS(i,j,bi,bj)  =
259         &           MAX(   R_low(i,j-1,bi,bj),   R_low(i,j,bi,bj) )
260             ENDDO
261            ENDDO
262           ENDDO
263          ENDDO
264          CALL EXCH_UV_XY_RS( rSurfW, rSurfS, .FALSE., myThid )
265          CALL EXCH_UV_XY_RS( rLowW,  rLowS,  .FALSE., myThid )
266    #endif /* ALLOW_NONHYDROSTATIC */
267    
268    C-    Write to disk: Total Column Thickness & hFac(C,W,S):
269    C     This I/O is now done in write_grid.F
270    c     CALL WRITE_FLD_XY_RS( 'Depth',' ',tmpfld,0,myThid)
271    c     CALL WRITE_FLD_XYZ_RS( 'hFacC',' ',hFacC,0,myThid)
272    c     CALL WRITE_FLD_XYZ_RS( 'hFacW',' ',hFacW,0,myThid)
273    c     CALL WRITE_FLD_XYZ_RS( 'hFacS',' ',hFacS,0,myThid)
274    
275          _BARRIER
276          CALL PLOT_FIELD_XYZRS( hFacC, 'hFacC' , Nr, 1, myThid )
277          CALL PLOT_FIELD_XYZRS( hFacW, 'hFacW' , Nr, 1, myThid )
278          CALL PLOT_FIELD_XYZRS( hFacS, 'hFacS' , Nr, 1, myThid )
279    
280  C     Masks and reciprocals of hFac[CWS]  C     Masks and reciprocals of hFac[CWS]
281        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
282         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
283          DO K=1,Nr          DO K=1,Nr
284           DO J=1,sNy           DO J=1-Oly,sNy+Oly
285            DO I=1,sNx            DO I=1-Olx,sNx+Olx
286             IF (HFacC(I,J,K,bi,bj) .NE. 0. _d 0 ) THEN             IF (hFacC(I,J,K,bi,bj) .NE. 0. ) THEN
287              recip_HFacC(I,J,K,bi,bj) = 1. _d 0 / HFacC(I,J,K,bi,bj)              recip_hFacC(I,J,K,bi,bj) = 1. _d 0 / hFacC(I,J,K,bi,bj)
288                maskC(I,J,K,bi,bj) = 1.
289             ELSE             ELSE
290              recip_HFacC(I,J,K,bi,bj) = 0. _d 0              recip_hFacC(I,J,K,bi,bj) = 0.
291                maskC(I,J,K,bi,bj) = 0.
292             ENDIF             ENDIF
293             IF (HFacW(I,J,K,bi,bj) .NE. 0. _d 0 ) THEN             IF (hFacW(I,J,K,bi,bj) .NE. 0. ) THEN
294              recip_HFacW(I,J,K,bi,bj) = 1. _d 0 / HFacW(I,J,K,bi,bj)              recip_hFacW(I,J,K,bi,bj) = 1. _d 0 / hFacW(I,J,K,bi,bj)
295              maskW(I,J,K,bi,bj) = 1. _d 0              maskW(I,J,K,bi,bj) = 1.
296             ELSE             ELSE
297              recip_HFacW(I,J,K,bi,bj) = 0. _d 0              recip_hFacW(I,J,K,bi,bj) = 0.
298              maskW(I,J,K,bi,bj) = 0.0 _d 0              maskW(I,J,K,bi,bj) = 0.
299             ENDIF             ENDIF
300             IF (HFacS(I,J,K,bi,bj) .NE. 0. _d 0 ) THEN             IF (hFacS(I,J,K,bi,bj) .NE. 0. ) THEN
301              recip_HFacS(I,J,K,bi,bj) = 1. _d 0 / HFacS(I,J,K,bi,bj)              recip_hFacS(I,J,K,bi,bj) = 1. _d 0 / hFacS(I,J,K,bi,bj)
302              maskS(I,J,K,bi,bj) = 1. _d 0              maskS(I,J,K,bi,bj) = 1.
303             ELSE             ELSE
304              recip_HFacS(I,J,K,bi,bj) = 0. _d 0              recip_hFacS(I,J,K,bi,bj) = 0.
305              maskS(I,J,K,bi,bj) = 0. _d 0              maskS(I,J,K,bi,bj) = 0.
306             ENDIF             ENDIF
307            ENDDO            ENDDO
308           ENDDO           ENDDO
309          ENDDO          ENDDO
310         ENDDO  C-    Calculate surface k index for interface W & S (U & V points)
311        ENDDO          DO J=1-Oly,sNy+Oly
312        _EXCH_XYZ_R4(recip_HFacC    , myThid )           DO I=1-Olx,sNx+Olx
313        _EXCH_XYZ_R4(recip_HFacW    , myThid )            ksurfW(I,J,bi,bj) = Nr+1
314        _EXCH_XYZ_R4(recip_HFacS    , myThid )            ksurfS(I,J,bi,bj) = Nr+1
315        _EXCH_XYZ_R4(maskW    , myThid )            DO k=Nr,1,-1
316        _EXCH_XYZ_R4(maskS    , myThid )             IF (hFacW(I,J,K,bi,bj).NE.0.) ksurfW(I,J,bi,bj) = k
317               IF (hFacS(I,J,K,bi,bj).NE.0.) ksurfS(I,J,bi,bj) = k
318  C     Calculate recipricols grid lengths            ENDDO
319        DO bj = myByLo(myThid), myByHi(myThid)            maskInW(I,J,bi,bj)= 0.
320         DO bi = myBxLo(myThid), myBxHi(myThid)            IF ( ksurfW(i,j,bi,bj).LE.Nr ) maskInW(I,J,bi,bj)= 1.
321          DO J=1,sNy            maskInS(I,J,bi,bj)= 0.
322           DO I=1,sNx            IF ( ksurfS(i,j,bi,bj).LE.Nr ) maskInS(I,J,bi,bj)= 1.
           recip_dxG(I,J,bi,bj)=1.d0/dxG(I,J,bi,bj)  
           recip_dyG(I,J,bi,bj)=1.d0/dyG(I,J,bi,bj)  
           recip_dxC(I,J,bi,bj)=1.d0/dxC(I,J,bi,bj)  
           recip_dyC(I,J,bi,bj)=1.d0/dyC(I,J,bi,bj)  
           recip_dxF(I,J,bi,bj)=1.d0/dxF(I,J,bi,bj)  
           recip_dyF(I,J,bi,bj)=1.d0/dyF(I,J,bi,bj)  
           recip_dxV(I,J,bi,bj)=1.d0/dxV(I,J,bi,bj)  
           recip_dyU(I,J,bi,bj)=1.d0/dyU(I,J,bi,bj)  
323           ENDDO           ENDDO
324          ENDDO          ENDDO
325    C - end bi,bj loops.
326         ENDDO         ENDDO
327        ENDDO        ENDDO
       _EXCH_XY_R4(recip_dxG, myThid )  
       _EXCH_XY_R4(recip_dyG, myThid )  
       _EXCH_XY_R4(recip_dxC, myThid )  
       _EXCH_XY_R4(recip_dyC, myThid )  
       _EXCH_XY_R4(recip_dxF, myThid )  
       _EXCH_XY_R4(recip_dyF, myThid )  
       _EXCH_XY_R4(recip_dxV, myThid )  
       _EXCH_XY_R4(recip_dyU, myThid )  
328    
329  #ifdef ALLOW_NONHYDROSTATIC  C     Calculate recipricols grid lengths
 C--   Calculate the reciprocal hfac distance/volume for W cells  
330        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
331         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
332          DO K=1,Nr          DO J=1-Oly,sNy+Oly
333           Km1=max(K-1,1)           DO I=1-Olx,sNx+Olx
334           hFacUpper=drF(Km1)/(drF(Km1)+drF(K))            IF ( dxG(I,J,bi,bj) .NE. 0. )
335           IF (Km1.EQ.K) hFacUpper=0.       &    recip_dxG(I,J,bi,bj)=1. _d 0/dxG(I,J,bi,bj)
336           hFacLower=drF(K)/(drF(Km1)+drF(K))            IF ( dyG(I,J,bi,bj) .NE. 0. )
337           DO J=1,sNy       &    recip_dyG(I,J,bi,bj)=1. _d 0/dyG(I,J,bi,bj)
338            DO I=1,sNx            IF ( dxC(I,J,bi,bj) .NE. 0. )
339             IF (hFacC(I,J,K,bi,bj).NE.0.) THEN       &    recip_dxC(I,J,bi,bj)=1. _d 0/dxC(I,J,bi,bj)
340              IF (hFacC(I,J,K,bi,bj).LE.0.5) THEN            IF ( dyC(I,J,bi,bj) .NE. 0. )
341              recip_hFacU(I,J,K,bi,bj)=       &    recip_dyC(I,J,bi,bj)=1. _d 0/dyC(I,J,bi,bj)
342       &         hFacUpper+hFacLower*hFacC(I,J,K,bi,bj)            IF ( dxF(I,J,bi,bj) .NE. 0. )
343              ELSE       &    recip_dxF(I,J,bi,bj)=1. _d 0/dxF(I,J,bi,bj)
344               recip_hFacU(I,J,K,bi,bj)=1.            IF ( dyF(I,J,bi,bj) .NE. 0. )
345              ENDIF       &    recip_dyF(I,J,bi,bj)=1. _d 0/dyF(I,J,bi,bj)
346             ELSE            IF ( dxV(I,J,bi,bj) .NE. 0. )
347              recip_hFacU(I,J,K,bi,bj)=0.       &    recip_dxV(I,J,bi,bj)=1. _d 0/dxV(I,J,bi,bj)
348             ENDIF            IF ( dyU(I,J,bi,bj) .NE. 0. )
349             IF (recip_hFacU(I,J,K,bi,bj).NE.0.)       &    recip_dyU(I,J,bi,bj)=1. _d 0/dyU(I,J,bi,bj)
350       &      recip_hFacU(I,J,K,bi,bj)=1./recip_hFacU(I,J,K,bi,bj)            IF ( rA(I,J,bi,bj) .NE. 0. )
351            ENDDO       &    recip_rA(I,J,bi,bj)=1. _d 0/rA(I,J,bi,bj)
352              IF ( rAs(I,J,bi,bj) .NE. 0. )
353         &    recip_rAs(I,J,bi,bj)=1. _d 0/rAs(I,J,bi,bj)
354              IF ( rAw(I,J,bi,bj) .NE. 0. )
355         &    recip_rAw(I,J,bi,bj)=1. _d 0/rAw(I,J,bi,bj)
356              IF ( rAz(I,J,bi,bj) .NE. 0. )
357         &    recip_rAz(I,J,bi,bj)=1. _d 0/rAz(I,J,bi,bj)
358           ENDDO           ENDDO
359          ENDDO          ENDDO
360         ENDDO         ENDDO
361        ENDDO        ENDDO
362        _EXCH_XY_R4(recip_hFacU, myThid )  
363  #endif  c #ifdef ALLOW_NONHYDROSTATIC
364  C  C--   Calculate "recip_hFacU" = reciprocal hfac distance/volume for W cells
365    C NOTE:  not used ; computed locally in CALC_GW
366    c #endif
367    
368        RETURN        RETURN
369        END        END
Legend:
Removed from v.1.17  
changed lines
  Added in v.1.44
  ViewVC Help
Powered by ViewVC 1.1.22