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