/[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.24 by cnh, Wed Sep 26 18:09:15 2001 UTC revision 1.37 by mlosch, Wed Sep 10 09:19:22 2008 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4    #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    
7  CBOP  CBOP
# Line 9  C     !INTERFACE: Line 10  C     !INTERFACE:
10        SUBROUTINE INI_MASKS_ETC( myThid )        SUBROUTINE INI_MASKS_ETC( myThid )
11  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
12  C     *==========================================================*  C     *==========================================================*
13  C     | SUBROUTINE INI_MASKS_ETC                                    C     | SUBROUTINE INI_MASKS_ETC
14  C     | o Initialise masks and topography factors                  C     | o Initialise masks and topography factors
15  C     *==========================================================*  C     *==========================================================*
16  C     | These arrays are used throughout the code and describe      C     | These arrays are used throughout the code and describe
17  C     | the topography of the domain through masks (0s and 1s)      C     | the topography of the domain through masks (0s and 1s)
18  C     | and fractional height factors (0<hFac<1). The latter        C     | and fractional height factors (0<hFac<1). The latter
19  C     | distinguish between the lopped-cell and full-step          C     | distinguish between the lopped-cell and full-step
20  C     | topographic representations.                                C     | topographic representations.
21  C     *==========================================================*  C     *==========================================================*
22  C     \ev  C     \ev
23    
# Line 28  C     === Global variables === Line 29  C     === Global variables ===
29  #include "PARAMS.h"  #include "PARAMS.h"
30  #include "GRID.h"  #include "GRID.h"
31  #include "SURFACE.h"  #include "SURFACE.h"
32    #ifdef ALLOW_SHELFICE
33    # include "SHELFICE.h"
34    #endif /* ALLOW_SHELFICE */
35    #ifdef ALLOW_EXCH2
36    # include "W2_EXCH2_TOPOLOGY.h"
37    # include "W2_EXCH2_PARAMS.h"
38    #endif
39    
40  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
41  C     == Routine arguments ==  C     == Routine arguments ==
42  C     myThid -  Number of this instance of INI_MASKS_ETC  C     myThid  ::  Number of this instance of INI_MASKS_ETC
43        INTEGER myThid        INTEGER myThid
44    
45  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
 C     == Local variables in common ==  
 C     tmpfld  - Temporary array used to compute & write Total Depth  
 C               has to be in common for multi threading  
       COMMON / LOCAL_INI_MASKS_ETC / tmpfld  
       _RS tmpfld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)  
46  C     == Local variables ==  C     == Local variables ==
47  C     bi,bj  - Loop counters  C     bi,bj   :: tile indices
48  C     I,J,K  C     I,J,K   :: Loop counters
49    C     tmpfld  :: Temporary array used to compute & write Total Depth
50          _RS tmpfld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
51        INTEGER bi, bj        INTEGER bi, bj
52        INTEGER  I, J, K        INTEGER I, J, K
 #ifdef ALLOW_NONHYDROSTATIC  
       INTEGER Km1  
       _RL hFacUpper,hFacLower  
 #endif  
53        _RL hFacCtmp        _RL hFacCtmp
54        _RL hFacMnSz        _RL hFacMnSz
55          _RL tileArea(nSx,nSy), threadArea
56    C     put tileArea in (local) common block to print from master-thread:
57          COMMON / LOCAL_INI_MASKS_ETC / tileArea
58          CHARACTER*(MAX_LEN_MBUF) msgBuf
59  CEOP  CEOP
60    
61  C- Calculate lopping factor hFacC : over-estimate the part inside of the domain  C- Calculate lopping factor hFacC : over-estimate the part inside of the domain
# Line 120  C      o Impose minimum fraction and/or Line 125  C      o Impose minimum fraction and/or
125            ENDDO            ENDDO
126           ENDDO           ENDDO
127          ENDDO          ENDDO
128           ENDDO
129          ENDDO
130    
131    #ifdef ALLOW_SHELFICE
132          IF ( useShelfIce ) THEN
133    C--   Modify lopping factor hFacC : Remove part outside of the domain
134    C     taking into account the Reference (=at rest) Surface Position Ro_shelfIce
135           CALL SHELFICE_UPDATE_MASKS(
136         I     rF, recip_drF,
137         U     hFacC,
138         I     myThid )
139          ENDIF
140    #endif /* ALLOW_SHELFICE */
141    
142  C-  Re-calculate Reference surface position, taking into account hFacC  C-  Re-calculate Reference surface position, taking into account hFacC
143  C   initialize Total column fluid thickness and surface k index  C   initialize Total column fluid thickness and surface k index
144  C       Note: if no fluid (continent) ==> ksurf = Nr+1  C       Note: if no fluid (continent) ==> ksurf = Nr+1
145          DO bj=myByLo(myThid), myByHi(myThid)
146           DO bi=myBxLo(myThid), myBxHi(myThid)
147          DO J=1-Oly,sNy+Oly          DO J=1-Oly,sNy+Oly
148           DO I=1-Olx,sNx+Olx           DO I=1-Olx,sNx+Olx
149            tmpfld(I,J,bi,bj) = 0.            tmpfld(I,J,bi,bj) = 0.
150            ksurfC(I,J,bi,bj) = Nr+1            ksurfC(I,J,bi,bj) = Nr+1
151              maskH(I,J,bi,bj) = 0.
152            Ro_surf(I,J,bi,bj) = R_low(I,J,bi,bj)            Ro_surf(I,J,bi,bj) = R_low(I,J,bi,bj)
153            DO K=Nr,1,-1            DO K=Nr,1,-1
154             Ro_surf(I,J,bi,bj) = Ro_surf(I,J,bi,bj)             Ro_surf(I,J,bi,bj) = Ro_surf(I,J,bi,bj)
155       &                        + drF(k)*hFacC(I,J,K,bi,bj)       &                        + drF(k)*hFacC(I,J,K,bi,bj)
156             IF (hFacC(I,J,K,bi,bj).NE.0.) THEN             IF (hFacC(I,J,K,bi,bj).NE.0.) THEN
157              ksurfC(I,J,bi,bj) = k              ksurfC(I,J,bi,bj) = k
158              tmpfld(i,j,bi,bj) = tmpfld(i,j,bi,bj) + 1.              maskH(I,J,bi,bj) = 1.
159                tmpfld(I,J,bi,bj) = tmpfld(I,J,bi,bj) + 1.
160               ENDIF
161              ENDDO
162              kLowC(I,J,bi,bj) = 0
163              DO K= 1, Nr
164               IF (hFacC(I,J,K,bi,bj).NE.0) THEN
165                  kLowC(I,J,bi,bj) = K
166             ENDIF             ENDIF
167            ENDDO            ENDDO
168           ENDDO           ENDDO
# Line 143  C - end bi,bj loops. Line 171  C - end bi,bj loops.
171         ENDDO         ENDDO
172        ENDDO        ENDDO
173    
174  C     CALL PLOT_FIELD_XYRS( tmpfld,  C     CALL PLOT_FIELD_XYRS( tmpfld,
175  C    &         'Model Depths K Index' , 1, myThid )  C    &         'Model Depths K Index' , 1, myThid )
176        CALL PLOT_FIELD_XYRS(R_low,        CALL PLOT_FIELD_XYRS(R_low,
177       &         'Model R_low (ini_masks_etc)', 1, myThid)       &         'Model R_low (ini_masks_etc)', 1, myThid)
178        CALL PLOT_FIELD_XYRS(Ro_surf,        CALL PLOT_FIELD_XYRS(Ro_surf,
179       &         'Model Ro_surf (ini_masks_etc)', 1, myThid)       &         'Model Ro_surf (ini_masks_etc)', 1, myThid)
180    
181  C     Calculate quantities derived from XY depth map  C     Calculate quantities derived from XY depth map
182          threadArea = 0. _d 0
183        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
184         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
185          DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
# Line 162  C         Inverse of fluid column thickn Line 191  C         Inverse of fluid column thickn
191            IF ( tmpfld(i,j,bi,bj) .LE. 0. ) THEN            IF ( tmpfld(i,j,bi,bj) .LE. 0. ) THEN
192             recip_Rcol(i,j,bi,bj) = 0.             recip_Rcol(i,j,bi,bj) = 0.
193            ELSE            ELSE
194             recip_Rcol(i,j,bi,bj) = 1. / tmpfld(i,j,bi,bj)             recip_Rcol(i,j,bi,bj) = 1. _d 0 / tmpfld(i,j,bi,bj)
195            ENDIF            ENDIF
196           ENDDO           ENDDO
197          ENDDO          ENDDO
198    C-      Compute the domain Area:
199            tileArea(bi,bj) = 0. _d 0
200            DO j=1,sNy
201             DO i=1,sNx
202              tileArea(bi,bj) = tileArea(bi,bj)
203         &                    + rA(i,j,bi,bj)*maskH(i,j,bi,bj)
204             ENDDO
205            ENDDO
206            threadArea = threadArea + tileArea(bi,bj)
207         ENDDO         ENDDO
208        ENDDO        ENDDO
209  C     _EXCH_XY_R4(   recip_Rcol, myThid )  C     _EXCH_XY_R4(   recip_Rcol, myThid )
210    #ifdef ALLOW_AUTODIFF_TAMC
211    C_jmc: apply GLOBAL_SUM to thread-local variable (not in common block)
212          _GLOBAL_SUM_R8( threadArea, myThid )
213    #else
214          CALL GLOBAL_SUM_TILE_RL( tileArea, threadArea, myThid )
215    #endif
216          _BEGIN_MASTER( myThid )
217          globalArea = threadArea
218    C-    list empty tiles:
219          msgBuf(1:1) = ' '
220          DO bj = 1,nSy
221           DO bi = 1,nSx
222            IF ( tileArea(bi,bj).EQ.0. _d 0 ) THEN
223    #ifdef ALLOW_EXCH2
224             WRITE(msgBuf,'(A,I6,A,I6,A)')
225         &     'Empty tile: #', W2_myTileList(bi), ' (bi=', bi,' )'
226    #else
227             WRITE(msgBuf,'(A,I6,I6)') 'Empty tile bi,bj=', bi, bj
228    #endif
229             CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
230         &                       SQUEEZE_RIGHT, myThid )
231            ENDIF
232           ENDDO
233          ENDDO
234          IF ( msgBuf(1:1).NE.' ' ) THEN
235             WRITE(msgBuf,'(A)') ' '
236             CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
237         &                       SQUEEZE_RIGHT, myThid )
238          ENDIF
239          _END_MASTER( myThid )
240    
241  C     hFacW and hFacS (at U and V points)  C     hFacW and hFacS (at U and V points)
242        DO bj=myByLo(myThid), myByHi(myThid)        DO bj=myByLo(myThid), myByHi(myThid)
243         DO bi=myBxLo(myThid), myBxHi(myThid)         DO bi=myBxLo(myThid), myBxHi(myThid)
244          DO K=1, Nr          DO K=1, Nr
245           DO J=1,sNy           DO J=1-Oly,sNy+Oly
246            DO I=1,sNx            hFacW(1-OLx,J,K,bi,bj)= 0.
247              DO I=2-Olx,sNx+Olx
248             hFacW(I,J,K,bi,bj)=             hFacW(I,J,K,bi,bj)=
249       &       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))
250              ENDDO
251             ENDDO
252             DO I=1-Olx,sNx+Olx
253               hFacS(I,1-OLy,K,bi,bj)= 0.
254             ENDDO
255             DO J=2-Oly,sNy+oly
256              DO I=1-Olx,sNx+Olx
257             hFacS(I,J,K,bi,bj)=             hFacS(I,J,K,bi,bj)=
258       &       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))
259            ENDDO            ENDDO
# Line 194  C     We should really supply a flag for Line 270  C     We should really supply a flag for
270          DO K=1, Nr          DO K=1, Nr
271           DO J=1-Oly,sNy+Oly           DO J=1-Oly,sNy+Oly
272            DO I=1-Olx,sNx+Olx            DO I=1-Olx,sNx+Olx
273             IF (DYG(I,J,bi,bj).EQ.0.) hFacW(I,J,K,bi,bj)=0.             IF (dyG(I,J,bi,bj).EQ.0.) hFacW(I,J,K,bi,bj)=0.
274             IF (DXG(I,J,bi,bj).EQ.0.) hFacS(I,J,K,bi,bj)=0.             IF (dxG(I,J,bi,bj).EQ.0.) hFacS(I,J,K,bi,bj)=0.
275            ENDDO            ENDDO
276           ENDDO           ENDDO
277          ENDDO          ENDDO
# Line 204  C     We should really supply a flag for Line 280  C     We should really supply a flag for
280    
281  C-    Write to disk: Total Column Thickness & hFac(C,W,S):  C-    Write to disk: Total Column Thickness & hFac(C,W,S):
282        _BARRIER        _BARRIER
283        _BEGIN_MASTER( myThid )  c     _BEGIN_MASTER( myThid )
284    C     This I/O is now done in write_grid.F
285  C     CALL MDSWRITEFIELD( 'Depth', writeBinaryPrec, .TRUE.,  C     CALL MDSWRITEFIELD( 'Depth', writeBinaryPrec, .TRUE.,
286  C    &                    'RS', 1, tmpfld, 1, -1, myThid )  C    &                    'RS', 1, tmpfld, 1, -1, myThid )
287        CALL WRITE_FLD_XY_RS( 'Depth',' ',tmpfld,0,myThid)  c     CALL WRITE_FLD_XY_RS( 'Depth',' ',tmpfld,0,myThid)
288        CALL WRITE_FLD_XYZ_RS( 'hFacC',' ',hFacC,0,myThid)  c     CALL WRITE_FLD_XYZ_RS( 'hFacC',' ',hFacC,0,myThid)
289        CALL WRITE_FLD_XYZ_RS( 'hFacW',' ',hFacW,0,myThid)  c     CALL WRITE_FLD_XYZ_RS( 'hFacW',' ',hFacW,0,myThid)
290        CALL WRITE_FLD_XYZ_RS( 'hFacS',' ',hFacS,0,myThid)  c     CALL WRITE_FLD_XYZ_RS( 'hFacS',' ',hFacS,0,myThid)
291        _END_MASTER(myThid)  c     _END_MASTER(myThid)
292    
293        CALL PLOT_FIELD_XYZRS( hFacC, 'hFacC' , Nr, 1, myThid )        CALL PLOT_FIELD_XYZRS( hFacC, 'hFacC' , Nr, 1, myThid )
294        CALL PLOT_FIELD_XYZRS( hFacW, 'hFacW' , Nr, 1, myThid )        CALL PLOT_FIELD_XYZRS( hFacW, 'hFacW' , Nr, 1, myThid )
# Line 223  C     Masks and reciprocals of hFac[CWS] Line 300  C     Masks and reciprocals of hFac[CWS]
300          DO K=1,Nr          DO K=1,Nr
301           DO J=1-Oly,sNy+Oly           DO J=1-Oly,sNy+Oly
302            DO I=1-Olx,sNx+Olx            DO I=1-Olx,sNx+Olx
303             IF (HFacC(I,J,K,bi,bj) .NE. 0. ) THEN             IF (hFacC(I,J,K,bi,bj) .NE. 0. ) THEN
304              recip_HFacC(I,J,K,bi,bj) = 1. / HFacC(I,J,K,bi,bj)              recip_hFacC(I,J,K,bi,bj) = 1. _d 0 / hFacC(I,J,K,bi,bj)
305              maskC(I,J,K,bi,bj) = 1.              maskC(I,J,K,bi,bj) = 1.
306             ELSE             ELSE
307              recip_HFacC(I,J,K,bi,bj) = 0.              recip_hFacC(I,J,K,bi,bj) = 0.
308              maskC(I,J,K,bi,bj) = 0.              maskC(I,J,K,bi,bj) = 0.
309             ENDIF             ENDIF
310             IF (HFacW(I,J,K,bi,bj) .NE. 0. ) THEN             IF (hFacW(I,J,K,bi,bj) .NE. 0. ) THEN
311              recip_HFacW(I,J,K,bi,bj) = 1. / HFacW(I,J,K,bi,bj)              recip_hFacW(I,J,K,bi,bj) = 1. _d 0 / hFacW(I,J,K,bi,bj)
312              maskW(I,J,K,bi,bj) = 1.              maskW(I,J,K,bi,bj) = 1.
313             ELSE             ELSE
314              recip_HFacW(I,J,K,bi,bj) = 0.              recip_hFacW(I,J,K,bi,bj) = 0.
315              maskW(I,J,K,bi,bj) = 0.              maskW(I,J,K,bi,bj) = 0.
316             ENDIF             ENDIF
317             IF (HFacS(I,J,K,bi,bj) .NE. 0. ) THEN             IF (hFacS(I,J,K,bi,bj) .NE. 0. ) THEN
318              recip_HFacS(I,J,K,bi,bj) = 1. / HFacS(I,J,K,bi,bj)              recip_hFacS(I,J,K,bi,bj) = 1. _d 0 / hFacS(I,J,K,bi,bj)
319              maskS(I,J,K,bi,bj) = 1.              maskS(I,J,K,bi,bj) = 1.
320             ELSE             ELSE
321              recip_HFacS(I,J,K,bi,bj) = 0.              recip_hFacS(I,J,K,bi,bj) = 0.
322              maskS(I,J,K,bi,bj) = 0.              maskS(I,J,K,bi,bj) = 0.
323             ENDIF             ENDIF
324            ENDDO            ENDDO
# Line 261  C-    Calculate surface k index for inte Line 338  C-    Calculate surface k index for inte
338  C - end bi,bj loops.  C - end bi,bj loops.
339         ENDDO         ENDDO
340        ENDDO        ENDDO
 C     _EXCH_XYZ_R4(recip_HFacC    , myThid )  
 C     _EXCH_XYZ_R4(recip_HFacW    , myThid )  
 C     _EXCH_XYZ_R4(recip_HFacS    , myThid )  
 C     _EXCH_XYZ_R4(maskW    , myThid )  
 C     _EXCH_XYZ_R4(maskS    , myThid )  
341    
342  C     Calculate recipricols grid lengths  C     Calculate recipricols grid lengths
343        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
# Line 273  C     Calculate recipricols grid lengths Line 345  C     Calculate recipricols grid lengths
345          DO J=1-Oly,sNy+Oly          DO J=1-Oly,sNy+Oly
346           DO I=1-Olx,sNx+Olx           DO I=1-Olx,sNx+Olx
347            IF ( dxG(I,J,bi,bj) .NE. 0. )            IF ( dxG(I,J,bi,bj) .NE. 0. )
348       &    recip_dxG(I,J,bi,bj)=1.d0/dxG(I,J,bi,bj)       &    recip_dxG(I,J,bi,bj)=1. _d 0/dxG(I,J,bi,bj)
349            IF ( dyG(I,J,bi,bj) .NE. 0. )            IF ( dyG(I,J,bi,bj) .NE. 0. )
350       &    recip_dyG(I,J,bi,bj)=1.d0/dyG(I,J,bi,bj)       &    recip_dyG(I,J,bi,bj)=1. _d 0/dyG(I,J,bi,bj)
351            IF ( dxC(I,J,bi,bj) .NE. 0. )            IF ( dxC(I,J,bi,bj) .NE. 0. )
352       &    recip_dxC(I,J,bi,bj)=1.d0/dxC(I,J,bi,bj)       &    recip_dxC(I,J,bi,bj)=1. _d 0/dxC(I,J,bi,bj)
353            IF ( dyC(I,J,bi,bj) .NE. 0. )            IF ( dyC(I,J,bi,bj) .NE. 0. )
354       &    recip_dyC(I,J,bi,bj)=1.d0/dyC(I,J,bi,bj)       &    recip_dyC(I,J,bi,bj)=1. _d 0/dyC(I,J,bi,bj)
355            IF ( dxF(I,J,bi,bj) .NE. 0. )            IF ( dxF(I,J,bi,bj) .NE. 0. )
356       &    recip_dxF(I,J,bi,bj)=1.d0/dxF(I,J,bi,bj)       &    recip_dxF(I,J,bi,bj)=1. _d 0/dxF(I,J,bi,bj)
357            IF ( dyF(I,J,bi,bj) .NE. 0. )            IF ( dyF(I,J,bi,bj) .NE. 0. )
358       &    recip_dyF(I,J,bi,bj)=1.d0/dyF(I,J,bi,bj)       &    recip_dyF(I,J,bi,bj)=1. _d 0/dyF(I,J,bi,bj)
359            IF ( dxV(I,J,bi,bj) .NE. 0. )            IF ( dxV(I,J,bi,bj) .NE. 0. )
360       &    recip_dxV(I,J,bi,bj)=1.d0/dxV(I,J,bi,bj)       &    recip_dxV(I,J,bi,bj)=1. _d 0/dxV(I,J,bi,bj)
361            IF ( dyU(I,J,bi,bj) .NE. 0. )            IF ( dyU(I,J,bi,bj) .NE. 0. )
362       &    recip_dyU(I,J,bi,bj)=1.d0/dyU(I,J,bi,bj)       &    recip_dyU(I,J,bi,bj)=1. _d 0/dyU(I,J,bi,bj)
363            IF ( rA(I,J,bi,bj) .NE. 0. )            IF ( rA(I,J,bi,bj) .NE. 0. )
364       &    recip_rA(I,J,bi,bj)=1.d0/rA(I,J,bi,bj)       &    recip_rA(I,J,bi,bj)=1. _d 0/rA(I,J,bi,bj)
365            IF ( rAs(I,J,bi,bj) .NE. 0. )            IF ( rAs(I,J,bi,bj) .NE. 0. )
366       &    recip_rAs(I,J,bi,bj)=1.d0/rAs(I,J,bi,bj)       &    recip_rAs(I,J,bi,bj)=1. _d 0/rAs(I,J,bi,bj)
367            IF ( rAw(I,J,bi,bj) .NE. 0. )            IF ( rAw(I,J,bi,bj) .NE. 0. )
368       &    recip_rAw(I,J,bi,bj)=1.d0/rAw(I,J,bi,bj)       &    recip_rAw(I,J,bi,bj)=1. _d 0/rAw(I,J,bi,bj)
369            IF ( rAz(I,J,bi,bj) .NE. 0. )            IF ( rAz(I,J,bi,bj) .NE. 0. )
370       &    recip_rAz(I,J,bi,bj)=1.d0/rAz(I,J,bi,bj)       &    recip_rAz(I,J,bi,bj)=1. _d 0/rAz(I,J,bi,bj)
371           ENDDO           ENDDO
372          ENDDO          ENDDO
373         ENDDO         ENDDO
374        ENDDO        ENDDO
 C     Do not need these since above denominators are valid over full range  
 C     _EXCH_XY_R4(recip_dxG, myThid )  
 C     _EXCH_XY_R4(recip_dyG, myThid )  
 C     _EXCH_XY_R4(recip_dxC, myThid )  
 C     _EXCH_XY_R4(recip_dyC, myThid )  
 C     _EXCH_XY_R4(recip_dxF, myThid )  
 C     _EXCH_XY_R4(recip_dyF, myThid )  
 C     _EXCH_XY_R4(recip_dxV, myThid )  
 C     _EXCH_XY_R4(recip_dyU, myThid )  
 C     _EXCH_XY_R4(recip_rAw, myThid )  
 C     _EXCH_XY_R4(recip_rAs, myThid )  
375    
376  #ifdef ALLOW_NONHYDROSTATIC  c #ifdef ALLOW_NONHYDROSTATIC
377  C--   Calculate the reciprocal hfac distance/volume for W cells  C--   Calculate "recip_hFacU" = reciprocal hfac distance/volume for W cells
378        DO bj = myByLo(myThid), myByHi(myThid)  C NOTE:  not used ; computed locally in CALC_GW
379         DO bi = myBxLo(myThid), myBxHi(myThid)  c #endif
380          DO K=1,Nr  
          Km1=max(K-1,1)  
          hFacUpper=drF(Km1)/(drF(Km1)+drF(K))  
          IF (Km1.EQ.K) hFacUpper=0.  
          hFacLower=drF(K)/(drF(Km1)+drF(K))  
          DO J=1-Oly,sNy+Oly  
           DO I=1-Olx,sNx+Olx  
            IF (hFacC(I,J,K,bi,bj).NE.0.) THEN  
             IF (hFacC(I,J,K,bi,bj).LE.0.5) THEN  
             recip_hFacU(I,J,K,bi,bj)=  
      &         hFacUpper+hFacLower*hFacC(I,J,K,bi,bj)  
             ELSE  
              recip_hFacU(I,J,K,bi,bj)=1.  
             ENDIF  
            ELSE  
             recip_hFacU(I,J,K,bi,bj)=0.  
            ENDIF  
            IF (recip_hFacU(I,J,K,bi,bj).NE.0.)  
      &      recip_hFacU(I,J,K,bi,bj)=1./recip_hFacU(I,J,K,bi,bj)  
           ENDDO  
          ENDDO  
         ENDDO  
        ENDDO  
       ENDDO  
 C     _EXCH_XY_R4(recip_hFacU, myThid )  
 #endif  
 C  
381        RETURN        RETURN
382        END        END
Legend:
Removed from v.1.24  
changed lines
  Added in v.1.37
  ViewVC Help
Powered by ViewVC 1.1.22