model/src/ini_masks_etc.F

C $Header: /u/gcmpack/models/MITgcmUV/model/src/ini_masks_etc.F,v 1.20.2.6 2001/05/01 13:21:40 adcroft Exp $
C $Name:  $

#include "CPP_OPTIONS.h"

CStartOfInterface
      SUBROUTINE INI_MASKS_ETC( myThid )
C     /==========================================================\
C     | SUBROUTINE INI_MASKS_ETC                                 |
C     | o Initialise masks and topography factors                |
C     |==========================================================|
C     | These arrays are used throughout the code and describe   |
C     | the topography of the domain through masks (0s and 1s)   |
C     | and fractional height factors (0<hFac<1). The latter     |
C     | distinguish between the lopped-cell and full-step        |
C     | topographic representations.                             |
C     \==========================================================/
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SURFACE.h"

C     == Routine arguments ==
C     myThid -  Number of this instance of INI_MASKS_ETC
      INTEGER myThid
CEndOfInterface

C     == Local variables ==
C     bi,bj  - Loop counters
C     I,J,K
      INTEGER bi, bj
      INTEGER  I, J, K
#ifdef ALLOW_NONHYDROSTATIC
      INTEGER Km1
      _RL hFacUpper,hFacLower
#endif
      _RL hFacCtmp
      _RL hFacMnSz
      _RS tmpfld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)

C- Calculate lopping factor hFacC : over-estimate the part inside of the domain
C    taking into account the lower_R Boundary (Bathymetrie / Top of Atmos)
      DO bj=myByLo(myThid), myByHi(myThid)
       DO bi=myBxLo(myThid), myBxHi(myThid)
        DO K=1, Nr
         hFacMnSz=max( hFacMin, min(hFacMinDr*recip_drF(k),1. _d 0) )
         DO J=1-Oly,sNy+Oly
          DO I=1-Olx,sNx+Olx
C      o Non-dimensional distance between grid bound. and domain lower_R bound.
           hFacCtmp = (rF(K)-R_low(I,J,bi,bj))*recip_drF(K)
C      o Select between, closed, open or partial (0,1,0-1)
           hFacCtmp=min( max( hFacCtmp, 0. _d 0) , 1. _d 0)
C      o Impose minimum fraction and/or size (dimensional)
           IF (hFacCtmp.LT.hFacMnSz) THEN
            IF (hFacCtmp.LT.hFacMnSz*0.5) THEN
             hFacC(I,J,K,bi,bj)=0.
            ELSE
             hFacC(I,J,K,bi,bj)=hFacMnSz
            ENDIF
           ELSE
             hFacC(I,J,K,bi,bj)=hFacCtmp
           ENDIF
          ENDDO
         ENDDO
        ENDDO

C-  Re-calculate lower-R Boundary position, taking into account hFacC
        DO J=1-Oly,sNy+Oly
         DO I=1-Olx,sNx+Olx
          R_low(I,J,bi,bj) = rF(1)
          DO K=Nr,1,-1
           R_low(I,J,bi,bj) = R_low(I,J,bi,bj)
     &                      - drF(k)*hFacC(I,J,K,bi,bj)
          ENDDO
         ENDDO
        ENDDO
C - end bi,bj loops.
       ENDDO
      ENDDO

C-  Calculate lopping factor hFacC : Remove part outside of the domain
C    taking into account the Reference (=at rest) Surface Position Ro_surf
      DO bj=myByLo(myThid), myByHi(myThid)
       DO bi=myBxLo(myThid), myBxHi(myThid)
        DO K=1, Nr
         hFacMnSz=max( hFacMin, min(hFacMinDr*recip_drF(k),1. _d 0) )
         DO J=1-Oly,sNy+Oly
          DO I=1-Olx,sNx+Olx
C      o Non-dimensional distance between grid boundary and model surface
           hFacCtmp = (rF(k)-Ro_surf(I,J,bi,bj))*recip_drF(K)
C      o Reduce the previous fraction : substract the outside part.
           hFacCtmp = hFacC(I,J,K,bi,bj) - max( hFacCtmp, 0. _d 0)
C      o set to zero if empty Column :
           hFacCtmp = max( hFacCtmp, 0. _d 0)
C      o Impose minimum fraction and/or size (dimensional)
           IF (hFacCtmp.LT.hFacMnSz) THEN
            IF (hFacCtmp.LT.hFacMnSz*0.5) THEN
             hFacC(I,J,K,bi,bj)=0.
            ELSE
             hFacC(I,J,K,bi,bj)=hFacMnSz
            ENDIF
           ELSE
             hFacC(I,J,K,bi,bj)=hFacCtmp
           ENDIF
          ENDDO
         ENDDO
        ENDDO

C-  Re-calculate Reference surface position, taking into account hFacC
C   initialize Total column fluid thickness and surface k index
        DO J=1-Oly,sNy+Oly
         DO I=1-Olx,sNx+Olx
          tmpfld(I,J,bi,bj) = 0.
          k_surf(I,J,bi,bj) = Nr
          Ro_surf(I,J,bi,bj) = R_low(I,J,bi,bj)
          DO K=Nr,1,-1
           Ro_surf(I,J,bi,bj) = Ro_surf(I,J,bi,bj)
     &                        + drF(k)*hFacC(I,J,K,bi,bj)
           IF (hFacC(I,J,K,bi,bj).NE.0.) THEN
            k_surf(I,J,bi,bj) = k
            tmpfld(i,j,bi,bj) = tmpfld(i,j,bi,bj) + 1.
           ENDIF
          ENDDO
         ENDDO
        ENDDO
C - end bi,bj loops.
       ENDDO
      ENDDO

C     CALL PLOT_FIELD_XYRS( tmpfld, 
C    &         'Model Depths K Index' , 1, myThid )
      CALL PLOT_FIELD_XYRS(R_low, 
     &         'Model R_low (ini_masks_etc)', 1, myThid)
      CALL PLOT_FIELD_XYRS(Ro_surf, 
     &         'Model Ro_surf (ini_masks_etc)', 1, myThid)

C     Calculate quantities derived from XY depth map
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO j=1-Oly,sNy+Oly
         DO i=1-Olx,sNx+Olx
C         Total fluid column thickness (r_unit) :
c         Rcolumn(i,j,bi,bj)= Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
          tmpfld(i,j,bi,bj) = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
C         Inverse of fluid column thickness (1/r_unit)
          IF ( tmpfld(i,j,bi,bj) .LE. 0. ) THEN
           recip_Rcol(i,j,bi,bj) = 0.
          ELSE
           recip_Rcol(i,j,bi,bj) = 1. / tmpfld(i,j,bi,bj)
          ENDIF
         ENDDO
        ENDDO
       ENDDO
      ENDDO
C     _EXCH_XY_R4(   recip_Rcol, myThid )
      CALL WRITE_FLD_XY_RS( 'Depth',' ',tmpfld,0,myThid)
      CALL WRITE_FLD_XYZ_RS( 'hFacC',' ',hFacC,0,myThid)
C     CALL MDSWRITEFIELD( 'Depth', writeBinaryPrec, .TRUE.,
C    &                    'RS', 1, tmpfld, 1, -1, myThid )

C     hFacW and hFacS (at U and V points)
      DO bj=myByLo(myThid), myByHi(myThid)
       DO bi=myBxLo(myThid), myBxHi(myThid)
        DO K=1, Nr
         DO J=1,sNy
          DO I=1,sNx
           hFacW(I,J,K,bi,bj)=
     &       MIN(hFacC(I,J,K,bi,bj),hFacC(I-1,J,K,bi,bj))
           hFacS(I,J,K,bi,bj)=
     &       MIN(hFacC(I,J,K,bi,bj),hFacC(I,J-1,K,bi,bj))
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      CALL EXCH_UV_XYZ_RS(hFacW,hFacS,.FALSE.,myThid)
C     The following block allows thin walls representation of non-periodic
C     boundaries such as happen on the lat-lon grid at the N/S poles.
C     We should really supply a flag for doing this.
      DO bj=myByLo(myThid), myByHi(myThid)
       DO bi=myBxLo(myThid), myBxHi(myThid)
        DO K=1, Nr
         DO J=1-Oly,sNy+Oly
          DO I=1-Olx,sNx+Olx
           IF (DYG(I,J,bi,bj).EQ.0.) hFacW(I,J,K,bi,bj)=0.
           IF (DXG(I,J,bi,bj).EQ.0.) hFacS(I,J,K,bi,bj)=0.
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO

      CALL PLOT_FIELD_XYZRS( hFacC, 'hFacC' , Nr, 1, myThid )
      CALL PLOT_FIELD_XYZRS( hFacW, 'hFacW' , Nr, 1, myThid )
      CALL PLOT_FIELD_XYZRS( hFacS, 'hFacS' , Nr, 1, myThid )

C     Masks and reciprocals of hFac[CWS]
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO K=1,Nr
         DO J=1-Oly,sNy+Oly
          DO I=1-Olx,sNx+Olx
           IF (HFacC(I,J,K,bi,bj) .NE. 0. ) THEN
            recip_HFacC(I,J,K,bi,bj) = 1. / HFacC(I,J,K,bi,bj)
            maskC(I,J,K,bi,bj) = 1.
           ELSE
            recip_HFacC(I,J,K,bi,bj) = 0.
            maskC(I,J,K,bi,bj) = 0.
           ENDIF
           IF (HFacW(I,J,K,bi,bj) .NE. 0. ) THEN
            recip_HFacW(I,J,K,bi,bj) = 1. / HFacW(I,J,K,bi,bj)
            maskW(I,J,K,bi,bj) = 1.
           ELSE
            recip_HFacW(I,J,K,bi,bj) = 0.
            maskW(I,J,K,bi,bj) = 0.
           ENDIF
           IF (HFacS(I,J,K,bi,bj) .NE. 0. ) THEN
            recip_HFacS(I,J,K,bi,bj) = 1. / HFacS(I,J,K,bi,bj)
            maskS(I,J,K,bi,bj) = 1.
           ELSE
            recip_HFacS(I,J,K,bi,bj) = 0.
            maskS(I,J,K,bi,bj) = 0.
           ENDIF
          ENDDO
         ENDDO
        ENDDO
       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 )

C     Calculate recipricols grid lengths
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO J=1-Oly,sNy+Oly
         DO I=1-Olx,sNx+Olx
          IF ( dxG(I,J,bi,bj) .NE. 0. )
     &    recip_dxG(I,J,bi,bj)=1.d0/dxG(I,J,bi,bj)
          IF ( dyG(I,J,bi,bj) .NE. 0. )
     &    recip_dyG(I,J,bi,bj)=1.d0/dyG(I,J,bi,bj)
          IF ( dxC(I,J,bi,bj) .NE. 0. )
     &    recip_dxC(I,J,bi,bj)=1.d0/dxC(I,J,bi,bj)
          IF ( dyC(I,J,bi,bj) .NE. 0. )
     &    recip_dyC(I,J,bi,bj)=1.d0/dyC(I,J,bi,bj)
          IF ( dxF(I,J,bi,bj) .NE. 0. )
     &    recip_dxF(I,J,bi,bj)=1.d0/dxF(I,J,bi,bj)
          IF ( dyF(I,J,bi,bj) .NE. 0. )
     &    recip_dyF(I,J,bi,bj)=1.d0/dyF(I,J,bi,bj)
          IF ( dxV(I,J,bi,bj) .NE. 0. )
     &    recip_dxV(I,J,bi,bj)=1.d0/dxV(I,J,bi,bj)
          IF ( dyU(I,J,bi,bj) .NE. 0. )
     &    recip_dyU(I,J,bi,bj)=1.d0/dyU(I,J,bi,bj)
          IF ( rA(I,J,bi,bj) .NE. 0. )
     &    recip_rA(I,J,bi,bj)=1.d0/rA(I,J,bi,bj)
          IF ( rAs(I,J,bi,bj) .NE. 0. )
     &    recip_rAs(I,J,bi,bj)=1.d0/rAs(I,J,bi,bj)
          IF ( rAw(I,J,bi,bj) .NE. 0. )
     &    recip_rAw(I,J,bi,bj)=1.d0/rAw(I,J,bi,bj)
          IF ( rAz(I,J,bi,bj) .NE. 0. )
     &    recip_rAz(I,J,bi,bj)=1.d0/rAz(I,J,bi,bj)
         ENDDO
        ENDDO
       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 )

#ifdef ALLOW_NONHYDROSTATIC
C--   Calculate the reciprocal hfac distance/volume for W cells
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        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
      RETURN
      END
1	adcroft	1.21	C $Header: /u/gcmpack/models/MITgcmUV/model/src/ini_masks_etc.F,v 1.20.2.6 2001/05/01 13:21:40 adcroft Exp $
2			C $Name: $
3	adcroft	1.1
4	cnh	1.11	#include "CPP_OPTIONS.h"
5	adcroft	1.1
6			CStartOfInterface
7			SUBROUTINE INI_MASKS_ETC( myThid )
8			C /==========================================================\
9			C \| SUBROUTINE INI_MASKS_ETC \|
10			C \| o Initialise masks and topography factors \|
11			C \|==========================================================\|
12			C \| These arrays are used throughout the code and describe \|
13			C \| the topography of the domain through masks (0s and 1s) \|
14			C \| and fractional height factors (0<hFac<1). The latter \|
15			C \| distinguish between the lopped-cell and full-step \|
16			C \| topographic representations. \|
17			C \==========================================================/
18	adcroft	1.13	IMPLICIT NONE
19	adcroft	1.1
20			C === Global variables ===
21			#include "SIZE.h"
22			#include "EEPARAMS.h"
23			#include "PARAMS.h"
24			#include "GRID.h"
25	adcroft	1.21	#include "SURFACE.h"
26	adcroft	1.1
27			C == Routine arguments ==
28	cnh	1.6	C myThid - Number of this instance of INI_MASKS_ETC
29	adcroft	1.1	INTEGER myThid
30			CEndOfInterface
31
32			C == Local variables ==
33			C bi,bj - Loop counters
34			C I,J,K
35			INTEGER bi, bj
36			INTEGER I, J, K
37	adcroft	1.15	#ifdef ALLOW_NONHYDROSTATIC
38			INTEGER Km1
39			_RL hFacUpper,hFacLower
40			#endif
41	adcroft	1.21	_RL hFacCtmp
42	adcroft	1.19	_RL hFacMnSz
43	adcroft	1.21	_RS tmpfld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
44	adcroft	1.19
45	adcroft	1.21	C- Calculate lopping factor hFacC : over-estimate the part inside of the domain
46			C taking into account the lower_R Boundary (Bathymetrie / Top of Atmos)
47	adcroft	1.2	DO bj=myByLo(myThid), myByHi(myThid)
48			DO bi=myBxLo(myThid), myBxHi(myThid)
49	cnh	1.4	DO K=1, Nr
50	adcroft	1.21	hFacMnSz=max( hFacMin, min(hFacMinDr*recip_drF(k),1. _d 0) )
51	adcroft	1.19	DO J=1-Oly,sNy+Oly
52			DO I=1-Olx,sNx+Olx
53	adcroft	1.21	C o Non-dimensional distance between grid bound. and domain lower_R bound.
54			hFacCtmp = (rF(K)-R_low(I,J,bi,bj))*recip_drF(K)
55			C o Select between, closed, open or partial (0,1,0-1)
56	adcroft	1.19	hFacCtmp=min( max( hFacCtmp, 0. _d 0) , 1. _d 0)
57	adcroft	1.21	C o Impose minimum fraction and/or size (dimensional)
58			IF (hFacCtmp.LT.hFacMnSz) THEN
59			IF (hFacCtmp.LT.hFacMnSz*0.5) THEN
60	adcroft	1.3	hFacC(I,J,K,bi,bj)=0.
61			ELSE
62	adcroft	1.19	hFacC(I,J,K,bi,bj)=hFacMnSz
63	adcroft	1.3	ENDIF
64	adcroft	1.21	ELSE
65			hFacC(I,J,K,bi,bj)=hFacCtmp
66	adcroft	1.2	ENDIF
67			ENDDO
68			ENDDO
69			ENDDO
70	adcroft	1.21
71			C- Re-calculate lower-R Boundary position, taking into account hFacC
72			DO J=1-Oly,sNy+Oly
73			DO I=1-Olx,sNx+Olx
74			R_low(I,J,bi,bj) = rF(1)
75			DO K=Nr,1,-1
76			R_low(I,J,bi,bj) = R_low(I,J,bi,bj)
77			& - drF(k)*hFacC(I,J,K,bi,bj)
78			ENDDO
79			ENDDO
80			ENDDO
81			C - end bi,bj loops.
82	adcroft	1.2	ENDDO
83			ENDDO
84	cnh	1.7
85	adcroft	1.21	C- Calculate lopping factor hFacC : Remove part outside of the domain
86			C taking into account the Reference (=at rest) Surface Position Ro_surf
87	adcroft	1.16	DO bj=myByLo(myThid), myByHi(myThid)
88			DO bi=myBxLo(myThid), myBxHi(myThid)
89	adcroft	1.21	DO K=1, Nr
90			hFacMnSz=max( hFacMin, min(hFacMinDr*recip_drF(k),1. _d 0) )
91			DO J=1-Oly,sNy+Oly
92			DO I=1-Olx,sNx+Olx
93			C o Non-dimensional distance between grid boundary and model surface
94			hFacCtmp = (rF(k)-Ro_surf(I,J,bi,bj))*recip_drF(K)
95			C o Reduce the previous fraction : substract the outside part.
96			hFacCtmp = hFacC(I,J,K,bi,bj) - max( hFacCtmp, 0. _d 0)
97			C o set to zero if empty Column :
98			hFacCtmp = max( hFacCtmp, 0. _d 0)
99			C o Impose minimum fraction and/or size (dimensional)
100			IF (hFacCtmp.LT.hFacMnSz) THEN
101			IF (hFacCtmp.LT.hFacMnSz*0.5) THEN
102			hFacC(I,J,K,bi,bj)=0.
103			ELSE
104			hFacC(I,J,K,bi,bj)=hFacMnSz
105			ENDIF
106			ELSE
107			hFacC(I,J,K,bi,bj)=hFacCtmp
108			ENDIF
109			ENDDO
110			ENDDO
111			ENDDO
112
113			C- Re-calculate Reference surface position, taking into account hFacC
114			C initialize Total column fluid thickness and surface k index
115	adcroft	1.19	DO J=1-Oly,sNy+Oly
116			DO I=1-Olx,sNx+Olx
117	adcroft	1.21	tmpfld(I,J,bi,bj) = 0.
118			k_surf(I,J,bi,bj) = Nr
119			Ro_surf(I,J,bi,bj) = R_low(I,J,bi,bj)
120			DO K=Nr,1,-1
121			Ro_surf(I,J,bi,bj) = Ro_surf(I,J,bi,bj)
122			& + drF(k)*hFacC(I,J,K,bi,bj)
123			IF (hFacC(I,J,K,bi,bj).NE.0.) THEN
124			k_surf(I,J,bi,bj) = k
125			tmpfld(i,j,bi,bj) = tmpfld(i,j,bi,bj) + 1.
126			ENDIF
127	adcroft	1.16	ENDDO
128			ENDDO
129			ENDDO
130	adcroft	1.21	C - end bi,bj loops.
131	adcroft	1.16	ENDDO
132			ENDDO
133	adcroft	1.21
134			C CALL PLOT_FIELD_XYRS( tmpfld,
135			C & 'Model Depths K Index' , 1, myThid )
136			CALL PLOT_FIELD_XYRS(R_low,
137			& 'Model R_low (ini_masks_etc)', 1, myThid)
138			CALL PLOT_FIELD_XYRS(Ro_surf,
139			& 'Model Ro_surf (ini_masks_etc)', 1, myThid)
140	adcroft	1.16
141			C Calculate quantities derived from XY depth map
142			DO bj = myByLo(myThid), myByHi(myThid)
143			DO bi = myBxLo(myThid), myBxHi(myThid)
144	adcroft	1.21	DO j=1-Oly,sNy+Oly
145			DO i=1-Olx,sNx+Olx
146			C Total fluid column thickness (r_unit) :
147			c Rcolumn(i,j,bi,bj)= Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
148			tmpfld(i,j,bi,bj) = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
149			C Inverse of fluid column thickness (1/r_unit)
150			IF ( tmpfld(i,j,bi,bj) .LE. 0. ) THEN
151			recip_Rcol(i,j,bi,bj) = 0.
152	adcroft	1.16	ELSE
153	adcroft	1.21	recip_Rcol(i,j,bi,bj) = 1. / tmpfld(i,j,bi,bj)
154	adcroft	1.16	ENDIF
155			ENDDO
156			ENDDO
157			ENDDO
158			ENDDO
159	adcroft	1.21	C _EXCH_XY_R4( recip_Rcol, myThid )
160			CALL WRITE_FLD_XY_RS( 'Depth',' ',tmpfld,0,myThid)
161			CALL WRITE_FLD_XYZ_RS( 'hFacC',' ',hFacC,0,myThid)
162			C CALL MDSWRITEFIELD( 'Depth', writeBinaryPrec, .TRUE.,
163			C & 'RS', 1, tmpfld, 1, -1, myThid )
164	adcroft	1.1
165			C hFacW and hFacS (at U and V points)
166			DO bj=myByLo(myThid), myByHi(myThid)
167			DO bi=myBxLo(myThid), myBxHi(myThid)
168	cnh	1.4	DO K=1, Nr
169	adcroft	1.1	DO J=1,sNy
170			DO I=1,sNx
171			hFacW(I,J,K,bi,bj)=
172			& MIN(hFacC(I,J,K,bi,bj),hFacC(I-1,J,K,bi,bj))
173			hFacS(I,J,K,bi,bj)=
174			& MIN(hFacC(I,J,K,bi,bj),hFacC(I,J-1,K,bi,bj))
175			ENDDO
176			ENDDO
177			ENDDO
178			ENDDO
179			ENDDO
180	adcroft	1.21	CALL EXCH_UV_XYZ_RS(hFacW,hFacS,.FALSE.,myThid)
181			C The following block allows thin walls representation of non-periodic
182			C boundaries such as happen on the lat-lon grid at the N/S poles.
183			C We should really supply a flag for doing this.
184	adcroft	1.19	DO bj=myByLo(myThid), myByHi(myThid)
185			DO bi=myBxLo(myThid), myBxHi(myThid)
186			DO K=1, Nr
187			DO J=1-Oly,sNy+Oly
188	adcroft	1.21	DO I=1-Olx,sNx+Olx
189			IF (DYG(I,J,bi,bj).EQ.0.) hFacW(I,J,K,bi,bj)=0.
190			IF (DXG(I,J,bi,bj).EQ.0.) hFacS(I,J,K,bi,bj)=0.
191	adcroft	1.19	ENDDO
192			ENDDO
193			ENDDO
194			ENDDO
195			ENDDO
196
197			CALL PLOT_FIELD_XYZRS( hFacC, 'hFacC' , Nr, 1, myThid )
198			CALL PLOT_FIELD_XYZRS( hFacW, 'hFacW' , Nr, 1, myThid )
199			CALL PLOT_FIELD_XYZRS( hFacS, 'hFacS' , Nr, 1, myThid )
200	adcroft	1.1
201			C Masks and reciprocals of hFac[CWS]
202			DO bj = myByLo(myThid), myByHi(myThid)
203			DO bi = myBxLo(myThid), myBxHi(myThid)
204	cnh	1.4	DO K=1,Nr
205	adcroft	1.19	DO J=1-Oly,sNy+Oly
206			DO I=1-Olx,sNx+Olx
207			IF (HFacC(I,J,K,bi,bj) .NE. 0. ) THEN
208			recip_HFacC(I,J,K,bi,bj) = 1. / HFacC(I,J,K,bi,bj)
209	adcroft	1.21	maskC(I,J,K,bi,bj) = 1.
210	adcroft	1.1	ELSE
211	adcroft	1.19	recip_HFacC(I,J,K,bi,bj) = 0.
212	adcroft	1.21	maskC(I,J,K,bi,bj) = 0.
213	adcroft	1.1	ENDIF
214	adcroft	1.19	IF (HFacW(I,J,K,bi,bj) .NE. 0. ) THEN
215			recip_HFacW(I,J,K,bi,bj) = 1. / HFacW(I,J,K,bi,bj)
216			maskW(I,J,K,bi,bj) = 1.
217	adcroft	1.1	ELSE
218	adcroft	1.19	recip_HFacW(I,J,K,bi,bj) = 0.
219			maskW(I,J,K,bi,bj) = 0.
220	adcroft	1.1	ENDIF
221	adcroft	1.19	IF (HFacS(I,J,K,bi,bj) .NE. 0. ) THEN
222			recip_HFacS(I,J,K,bi,bj) = 1. / HFacS(I,J,K,bi,bj)
223			maskS(I,J,K,bi,bj) = 1.
224	adcroft	1.1	ELSE
225	adcroft	1.19	recip_HFacS(I,J,K,bi,bj) = 0.
226			maskS(I,J,K,bi,bj) = 0.
227	adcroft	1.1	ENDIF
228			ENDDO
229			ENDDO
230			ENDDO
231			ENDDO
232			ENDDO
233	adcroft	1.19	C _EXCH_XYZ_R4(recip_HFacC , myThid )
234			C _EXCH_XYZ_R4(recip_HFacW , myThid )
235			C _EXCH_XYZ_R4(recip_HFacS , myThid )
236			C _EXCH_XYZ_R4(maskW , myThid )
237			C _EXCH_XYZ_R4(maskS , myThid )
238	adcroft	1.1
239			C Calculate recipricols grid lengths
240			DO bj = myByLo(myThid), myByHi(myThid)
241			DO bi = myBxLo(myThid), myBxHi(myThid)
242	adcroft	1.19	DO J=1-Oly,sNy+Oly
243			DO I=1-Olx,sNx+Olx
244			IF ( dxG(I,J,bi,bj) .NE. 0. )
245			& recip_dxG(I,J,bi,bj)=1.d0/dxG(I,J,bi,bj)
246			IF ( dyG(I,J,bi,bj) .NE. 0. )
247			& recip_dyG(I,J,bi,bj)=1.d0/dyG(I,J,bi,bj)
248			IF ( dxC(I,J,bi,bj) .NE. 0. )
249			& recip_dxC(I,J,bi,bj)=1.d0/dxC(I,J,bi,bj)
250			IF ( dyC(I,J,bi,bj) .NE. 0. )
251			& recip_dyC(I,J,bi,bj)=1.d0/dyC(I,J,bi,bj)
252			IF ( dxF(I,J,bi,bj) .NE. 0. )
253			& recip_dxF(I,J,bi,bj)=1.d0/dxF(I,J,bi,bj)
254			IF ( dyF(I,J,bi,bj) .NE. 0. )
255			& recip_dyF(I,J,bi,bj)=1.d0/dyF(I,J,bi,bj)
256			IF ( dxV(I,J,bi,bj) .NE. 0. )
257			& recip_dxV(I,J,bi,bj)=1.d0/dxV(I,J,bi,bj)
258			IF ( dyU(I,J,bi,bj) .NE. 0. )
259			& recip_dyU(I,J,bi,bj)=1.d0/dyU(I,J,bi,bj)
260			IF ( rA(I,J,bi,bj) .NE. 0. )
261			& recip_rA(I,J,bi,bj)=1.d0/rA(I,J,bi,bj)
262			IF ( rAs(I,J,bi,bj) .NE. 0. )
263			& recip_rAs(I,J,bi,bj)=1.d0/rAs(I,J,bi,bj)
264			IF ( rAw(I,J,bi,bj) .NE. 0. )
265			& recip_rAw(I,J,bi,bj)=1.d0/rAw(I,J,bi,bj)
266			IF ( rAz(I,J,bi,bj) .NE. 0. )
267			& recip_rAz(I,J,bi,bj)=1.d0/rAz(I,J,bi,bj)
268	adcroft	1.1	ENDDO
269			ENDDO
270			ENDDO
271			ENDDO
272	adcroft	1.19	C Do not need these since above denominators are valid over full range
273			C _EXCH_XY_R4(recip_dxG, myThid )
274			C _EXCH_XY_R4(recip_dyG, myThid )
275			C _EXCH_XY_R4(recip_dxC, myThid )
276			C _EXCH_XY_R4(recip_dyC, myThid )
277			C _EXCH_XY_R4(recip_dxF, myThid )
278			C _EXCH_XY_R4(recip_dyF, myThid )
279			C _EXCH_XY_R4(recip_dxV, myThid )
280			C _EXCH_XY_R4(recip_dyU, myThid )
281			C _EXCH_XY_R4(recip_rAw, myThid )
282			C _EXCH_XY_R4(recip_rAs, myThid )
283	adcroft	1.1
284	adcroft	1.15	#ifdef ALLOW_NONHYDROSTATIC
285			C-- Calculate the reciprocal hfac distance/volume for W cells
286			DO bj = myByLo(myThid), myByHi(myThid)
287			DO bi = myBxLo(myThid), myBxHi(myThid)
288			DO K=1,Nr
289			Km1=max(K-1,1)
290			hFacUpper=drF(Km1)/(drF(Km1)+drF(K))
291			IF (Km1.EQ.K) hFacUpper=0.
292			hFacLower=drF(K)/(drF(Km1)+drF(K))
293	adcroft	1.19	DO J=1-Oly,sNy+Oly
294			DO I=1-Olx,sNx+Olx
295	adcroft	1.15	IF (hFacC(I,J,K,bi,bj).NE.0.) THEN
296			IF (hFacC(I,J,K,bi,bj).LE.0.5) THEN
297			recip_hFacU(I,J,K,bi,bj)=
298			& hFacUpper+hFacLower*hFacC(I,J,K,bi,bj)
299			ELSE
300			recip_hFacU(I,J,K,bi,bj)=1.
301			ENDIF
302			ELSE
303			recip_hFacU(I,J,K,bi,bj)=0.
304			ENDIF
305			IF (recip_hFacU(I,J,K,bi,bj).NE.0.)
306			& recip_hFacU(I,J,K,bi,bj)=1./recip_hFacU(I,J,K,bi,bj)
307			ENDDO
308			ENDDO
309			ENDDO
310			ENDDO
311			ENDDO
312	adcroft	1.19	C _EXCH_XY_R4(recip_hFacU, myThid )
313	adcroft	1.15	#endif
314	adcroft	1.1	C
315			RETURN
316			END