model/src/ini_masks_etc.F

C $Header: /u/gcmpack/MITgcm/model/src/ini_masks_etc.F,v 1.31 2006/02/07 11:47:48 mlosch Exp $
C $Name:  $

#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: INI_MASKS_ETC
C     !INTERFACE:
      SUBROUTINE INI_MASKS_ETC( myThid )
C     !DESCRIPTION: \bv
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     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SURFACE.h"
#ifdef ALLOW_SHELFICE
# include "SHELFICE.h"
#endif /* ALLOW_SHELFICE */

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     myThid  ::  Number of this instance of INI_MASKS_ETC
      INTEGER myThid

C     !LOCAL VARIABLES:
C     == Local variables ==
C     bi,bj   :: tile indices
C     I,J,K   :: Loop counters
C     tmpfld  :: Temporary array used to compute & write Total Depth
      _RS tmpfld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      INTEGER bi, bj
      INTEGER I, J, K
      _RL hFacCtmp
      _RL hFacMnSz
      _RL tileArea
CEOP

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

#ifdef ALLOW_SHELFICE
C--   compute contributions of shelf ice to looping factors
        IF ( useShelfIce ) THEN
        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)-R_shelfIce(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
        ENDIF
#endif /* ALLOW_SHELFICE */

C-  Re-calculate Reference surface position, taking into account hFacC
C   initialize Total column fluid thickness and surface k index
C       Note: if no fluid (continent) ==> ksurf = Nr+1
        DO J=1-Oly,sNy+Oly
         DO I=1-Olx,sNx+Olx
          tmpfld(I,J,bi,bj) = 0.
          ksurfC(I,J,bi,bj) = Nr+1
          maskH(I,J,bi,bj) = 0.
          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
            ksurfC(I,J,bi,bj) = k
            maskH(I,J,bi,bj) = 1.
            tmpfld(I,J,bi,bj) = tmpfld(I,J,bi,bj) + 1.
           ENDIF
          ENDDO
          kLowC(I,J,bi,bj) = 0
          DO K= 1, Nr
           IF (hFacC(I,J,K,bi,bj).NE.0) THEN
              kLowC(I,J,bi,bj) = K
           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
      globalArea = 0. _d 0
      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. _d 0 / tmpfld(i,j,bi,bj)
          ENDIF
         ENDDO
        ENDDO
C-      Compute the domain Area:
        tileArea = 0. _d 0
        DO j=1,sNy
         DO i=1,sNx
          tileArea = tileArea + rA(i,j,bi,bj)*maskH(i,j,bi,bj)
         ENDDO
        ENDDO
        globalArea = globalArea + tileArea
       ENDDO
      ENDDO
C     _EXCH_XY_R4(   recip_Rcol, myThid )
      _GLOBAL_SUM_R8( globalArea, 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-Oly,sNy+Oly
          DO I=2-Olx,sNx+Olx
           hFacW(I,J,K,bi,bj)=
     &       MIN(hFacC(I,J,K,bi,bj),hFacC(I-1,J,K,bi,bj))
          ENDDO
         ENDDO
         DO J=2-Oly,sNy+oly
          DO I=1-Olx,sNx+Olx
           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

C-    Write to disk: Total Column Thickness & hFac(C,W,S):
      _BARRIER
c     _BEGIN_MASTER( myThid )
C     This I/O is now done in write_grid.F
C     CALL MDSWRITEFIELD( 'Depth', writeBinaryPrec, .TRUE.,
C    &                    'RS', 1, tmpfld, 1, -1, myThid )
c     CALL WRITE_FLD_XY_RS( 'Depth',' ',tmpfld,0,myThid)
c     CALL WRITE_FLD_XYZ_RS( 'hFacC',' ',hFacC,0,myThid)
c     CALL WRITE_FLD_XYZ_RS( 'hFacW',' ',hFacW,0,myThid)
c     CALL WRITE_FLD_XYZ_RS( 'hFacS',' ',hFacS,0,myThid)
c     _END_MASTER(myThid)

      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. _d 0 / 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. _d 0 / 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. _d 0 / 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
C-    Calculate surface k index for interface W & S (U & V points)
        DO J=1-Oly,sNy+Oly
         DO I=1-Olx,sNx+Olx
          ksurfW(I,J,bi,bj) = Nr+1
          ksurfS(I,J,bi,bj) = Nr+1
          DO k=Nr,1,-1
           IF (hFacW(I,J,K,bi,bj).NE.0.) ksurfW(I,J,bi,bj) = k
           IF (hFacS(I,J,K,bi,bj).NE.0.) ksurfS(I,J,bi,bj) = k
          ENDDO
         ENDDO
        ENDDO
C - end bi,bj loops.
       ENDDO
      ENDDO

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. _d 0/dxG(I,J,bi,bj)
          IF ( dyG(I,J,bi,bj) .NE. 0. )
     &    recip_dyG(I,J,bi,bj)=1. _d 0/dyG(I,J,bi,bj)
          IF ( dxC(I,J,bi,bj) .NE. 0. )
     &    recip_dxC(I,J,bi,bj)=1. _d 0/dxC(I,J,bi,bj)
          IF ( dyC(I,J,bi,bj) .NE. 0. )
     &    recip_dyC(I,J,bi,bj)=1. _d 0/dyC(I,J,bi,bj)
          IF ( dxF(I,J,bi,bj) .NE. 0. )
     &    recip_dxF(I,J,bi,bj)=1. _d 0/dxF(I,J,bi,bj)
          IF ( dyF(I,J,bi,bj) .NE. 0. )
     &    recip_dyF(I,J,bi,bj)=1. _d 0/dyF(I,J,bi,bj)
          IF ( dxV(I,J,bi,bj) .NE. 0. )
     &    recip_dxV(I,J,bi,bj)=1. _d 0/dxV(I,J,bi,bj)
          IF ( dyU(I,J,bi,bj) .NE. 0. )
     &    recip_dyU(I,J,bi,bj)=1. _d 0/dyU(I,J,bi,bj)
          IF ( rA(I,J,bi,bj) .NE. 0. )
     &    recip_rA(I,J,bi,bj)=1. _d 0/rA(I,J,bi,bj)
          IF ( rAs(I,J,bi,bj) .NE. 0. )
     &    recip_rAs(I,J,bi,bj)=1. _d 0/rAs(I,J,bi,bj)
          IF ( rAw(I,J,bi,bj) .NE. 0. )
     &    recip_rAw(I,J,bi,bj)=1. _d 0/rAw(I,J,bi,bj)
          IF ( rAz(I,J,bi,bj) .NE. 0. )
     &    recip_rAz(I,J,bi,bj)=1. _d 0/rAz(I,J,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO

c #ifdef ALLOW_NONHYDROSTATIC
C--   Calculate "recip_hFacU" = reciprocal hfac distance/volume for W cells
C     not used ; computed locally in CALC_GW
c #endif

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