model/src/ini_masks_etc.F

C $Header: /u/gcmpack/MITgcm/model/src/ini_masks_etc.F,v 1.39 2009/02/27 01:05:23 jmc 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_EXCH2
# include "W2_EXCH2_TOPOLOGY.h"
# include "W2_EXCH2_PARAMS.h"
#endif /* ALLOW_EXCH2 */

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(nSx,nSy), threadArea
C     put tileArea in (local) common block to print from master-thread:
      COMMON / LOCAL_INI_MASKS_ETC / tileArea
      CHARACTER*(MAX_LEN_MBUF) msgBuf
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
       ENDDO
      ENDDO

#ifdef ALLOW_SHELFICE
      IF ( useShelfIce ) THEN
C--   Modify lopping factor hFacC : Remove part outside of the domain
C     taking into account the Reference (=at rest) Surface Position Ro_shelfIce
       CALL SHELFICE_UPDATE_MASKS(
     I     rF, recip_drF,
     U     hFacC,
     I     myThid )
      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 bj=myByLo(myThid), myByHi(myThid)
       DO bi=myBxLo(myThid), myBxHi(myThid)
        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
      threadArea = 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(bi,bj) = 0. _d 0
        DO j=1,sNy
         DO i=1,sNx
          tileArea(bi,bj) = tileArea(bi,bj)
     &                    + rA(i,j,bi,bj)*maskH(i,j,bi,bj)
         ENDDO
        ENDDO
        threadArea = threadArea + tileArea(bi,bj)
       ENDDO
      ENDDO
C     _EXCH_XY_RS(   recip_Rcol, myThid )
#ifdef ALLOW_AUTODIFF_TAMC
C_jmc: apply GLOBAL_SUM to thread-local variable (not in common block)
      _GLOBAL_SUM_RL( threadArea, myThid )
#else
      CALL GLOBAL_SUM_TILE_RL( tileArea, threadArea, myThid )
#endif
      _BEGIN_MASTER( myThid )
      globalArea = threadArea
C-    list empty tiles:
      msgBuf(1:1) = ' '
      DO bj = 1,nSy
       DO bi = 1,nSx
        IF ( tileArea(bi,bj).EQ.0. _d 0 ) THEN
#ifdef ALLOW_EXCH2
         WRITE(msgBuf,'(A,I6,A,I6,A)')
     &     'Empty tile: #', W2_myTileList(bi), ' (bi=', bi,' )'
#else
         WRITE(msgBuf,'(A,I6,I6)') 'Empty tile bi,bj=', bi, bj
#endif
         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                       SQUEEZE_RIGHT, myThid )
        ENDIF
       ENDDO
      ENDDO
      IF ( msgBuf(1:1).NE.' ' ) THEN
         WRITE(msgBuf,'(A)') ' '
         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                       SQUEEZE_RIGHT, myThid )
      ENDIF
      _END_MASTER( 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
          hFacW(1-OLx,J,K,bi,bj)= 0.
          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 I=1-Olx,sNx+Olx
           hFacS(I,1-OLy,K,bi,bj)= 0.
         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

#ifdef ALLOW_NONHYDROSTATIC
C     rLow & reference rSurf at Western & Southern edges (U and V points)
      DO bj=myByLo(myThid), myByHi(myThid)
       DO bi=myBxLo(myThid), myBxHi(myThid)
        I = 1-OlX
        DO J=1-Oly,sNy+Oly
           rLowW (i,j,bi,bj) = 0.
           rSurfW(i,j,bi,bj) = 0.
        ENDDO
        J = 1-Oly
        DO I=1-Olx,sNx+Olx
           rLowS (i,j,bi,bj) = 0.
           rSurfS(i,j,bi,bj) = 0.
        ENDDO
        DO J=1-Oly,sNy+Oly
         DO I=2-Olx,sNx+Olx
           rSurfW(i,j,bi,bj) =
     &           MIN( Ro_surf(i-1,j,bi,bj), Ro_surf(i,j,bi,bj) )
           rLowW(i,j,bi,bj)  =
     &           MAX(   R_low(i-1,j,bi,bj),   R_low(i,j,bi,bj) )
         ENDDO
        ENDDO
        DO J=2-Oly,sNy+Oly
         DO I=1-Olx,sNx+Olx
           rSurfS(i,j,bi,bj) =
     &           MIN( Ro_surf(i,j-1,bi,bj), Ro_surf(i,j,bi,bj) )
           rLowS(i,j,bi,bj)  =
     &           MAX(   R_low(i,j-1,bi,bj),   R_low(i,j,bi,bj) )
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      CALL EXCH_UV_XY_RS( rSurfW, rSurfS, .FALSE., myThid )
      CALL EXCH_UV_XY_RS( rLowW,  rLowS,  .FALSE., myThid )
#endif /* ALLOW_NONHYDROSTATIC */

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 NOTE:  not used ; computed locally in CALC_GW
c #endif

      RETURN
      END
1	jmc	1.40	C $Header: /u/gcmpack/MITgcm/model/src/ini_masks_etc.F,v 1.39 2009/02/27 01:05:23 jmc 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	jmc	1.36	#ifdef ALLOW_EXCH2
33			# include "W2_EXCH2_TOPOLOGY.h"
34			# include "W2_EXCH2_PARAMS.h"
35	mlosch	1.38	#endif /* ALLOW_EXCH2 */
36	adcroft	1.1
37	cnh	1.24	C !INPUT/OUTPUT PARAMETERS:
38	adcroft	1.1	C == Routine arguments ==
39	jmc	1.32	C myThid :: Number of this instance of INI_MASKS_ETC
40	adcroft	1.1	INTEGER myThid
41
42	cnh	1.24	C !LOCAL VARIABLES:
43	jmc	1.32	C == Local variables ==
44			C bi,bj :: tile indices
45			C I,J,K :: Loop counters
46			C tmpfld :: Temporary array used to compute & write Total Depth
47	jmc	1.22	_RS tmpfld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
48	adcroft	1.1	INTEGER bi, bj
49	jmc	1.32	INTEGER I, J, K
50	adcroft	1.21	_RL hFacCtmp
51	adcroft	1.19	_RL hFacMnSz
52	jmc	1.35	_RL tileArea(nSx,nSy), threadArea
53	jmc	1.36	C put tileArea in (local) common block to print from master-thread:
54			COMMON / LOCAL_INI_MASKS_ETC / tileArea
55			CHARACTER*(MAX_LEN_MBUF) msgBuf
56	cnh	1.24	CEOP
57	adcroft	1.19
58	adcroft	1.21	C- Calculate lopping factor hFacC : over-estimate the part inside of the domain
59			C taking into account the lower_R Boundary (Bathymetrie / Top of Atmos)
60	adcroft	1.2	DO bj=myByLo(myThid), myByHi(myThid)
61			DO bi=myBxLo(myThid), myBxHi(myThid)
62	cnh	1.4	DO K=1, Nr
63	adcroft	1.21	hFacMnSz=max( hFacMin, min(hFacMinDr*recip_drF(k),1. _d 0) )
64	adcroft	1.19	DO J=1-Oly,sNy+Oly
65			DO I=1-Olx,sNx+Olx
66	adcroft	1.21	C o Non-dimensional distance between grid bound. and domain lower_R bound.
67			hFacCtmp = (rF(K)-R_low(I,J,bi,bj))*recip_drF(K)
68			C o Select between, closed, open or partial (0,1,0-1)
69	adcroft	1.19	hFacCtmp=min( max( hFacCtmp, 0. _d 0) , 1. _d 0)
70	adcroft	1.21	C o Impose minimum fraction and/or size (dimensional)
71			IF (hFacCtmp.LT.hFacMnSz) THEN
72			IF (hFacCtmp.LT.hFacMnSz*0.5) THEN
73	adcroft	1.3	hFacC(I,J,K,bi,bj)=0.
74			ELSE
75	adcroft	1.19	hFacC(I,J,K,bi,bj)=hFacMnSz
76	adcroft	1.3	ENDIF
77	adcroft	1.21	ELSE
78			hFacC(I,J,K,bi,bj)=hFacCtmp
79	adcroft	1.2	ENDIF
80			ENDDO
81			ENDDO
82			ENDDO
83	adcroft	1.21
84			C- Re-calculate lower-R Boundary position, taking into account hFacC
85			DO J=1-Oly,sNy+Oly
86			DO I=1-Olx,sNx+Olx
87			R_low(I,J,bi,bj) = rF(1)
88			DO K=Nr,1,-1
89			R_low(I,J,bi,bj) = R_low(I,J,bi,bj)
90			& - drF(k)*hFacC(I,J,K,bi,bj)
91			ENDDO
92			ENDDO
93			ENDDO
94			C - end bi,bj loops.
95	adcroft	1.2	ENDDO
96			ENDDO
97	cnh	1.7
98	adcroft	1.21	C- Calculate lopping factor hFacC : Remove part outside of the domain
99			C taking into account the Reference (=at rest) Surface Position Ro_surf
100	adcroft	1.16	DO bj=myByLo(myThid), myByHi(myThid)
101			DO bi=myBxLo(myThid), myBxHi(myThid)
102	adcroft	1.21	DO K=1, Nr
103			hFacMnSz=max( hFacMin, min(hFacMinDr*recip_drF(k),1. _d 0) )
104			DO J=1-Oly,sNy+Oly
105			DO I=1-Olx,sNx+Olx
106			C o Non-dimensional distance between grid boundary and model surface
107			hFacCtmp = (rF(k)-Ro_surf(I,J,bi,bj))*recip_drF(K)
108			C o Reduce the previous fraction : substract the outside part.
109			hFacCtmp = hFacC(I,J,K,bi,bj) - max( hFacCtmp, 0. _d 0)
110			C o set to zero if empty Column :
111			hFacCtmp = max( hFacCtmp, 0. _d 0)
112			C o Impose minimum fraction and/or size (dimensional)
113			IF (hFacCtmp.LT.hFacMnSz) THEN
114			IF (hFacCtmp.LT.hFacMnSz*0.5) THEN
115			hFacC(I,J,K,bi,bj)=0.
116			ELSE
117			hFacC(I,J,K,bi,bj)=hFacMnSz
118			ENDIF
119			ELSE
120			hFacC(I,J,K,bi,bj)=hFacCtmp
121			ENDIF
122			ENDDO
123			ENDDO
124			ENDDO
125	mlosch	1.37	ENDDO
126			ENDDO
127	adcroft	1.21
128	mlosch	1.31	#ifdef ALLOW_SHELFICE
129	mlosch	1.37	IF ( useShelfIce ) THEN
130			C-- Modify lopping factor hFacC : Remove part outside of the domain
131			C taking into account the Reference (=at rest) Surface Position Ro_shelfIce
132	jmc	1.39	CALL SHELFICE_UPDATE_MASKS(
133	mlosch	1.37	I rF, recip_drF,
134	jmc	1.39	U hFacC,
135	mlosch	1.37	I myThid )
136			ENDIF
137	mlosch	1.31	#endif /* ALLOW_SHELFICE */
138
139	adcroft	1.21	C- Re-calculate Reference surface position, taking into account hFacC
140			C initialize Total column fluid thickness and surface k index
141	jmc	1.23	C Note: if no fluid (continent) ==> ksurf = Nr+1
142	mlosch	1.37	DO bj=myByLo(myThid), myByHi(myThid)
143			DO bi=myBxLo(myThid), myBxHi(myThid)
144	adcroft	1.19	DO J=1-Oly,sNy+Oly
145			DO I=1-Olx,sNx+Olx
146	adcroft	1.21	tmpfld(I,J,bi,bj) = 0.
147	jmc	1.23	ksurfC(I,J,bi,bj) = Nr+1
148	jmc	1.32	maskH(I,J,bi,bj) = 0.
149	adcroft	1.21	Ro_surf(I,J,bi,bj) = R_low(I,J,bi,bj)
150			DO K=Nr,1,-1
151			Ro_surf(I,J,bi,bj) = Ro_surf(I,J,bi,bj)
152			& + drF(k)*hFacC(I,J,K,bi,bj)
153			IF (hFacC(I,J,K,bi,bj).NE.0.) THEN
154	jmc	1.23	ksurfC(I,J,bi,bj) = k
155	jmc	1.32	maskH(I,J,bi,bj) = 1.
156			tmpfld(I,J,bi,bj) = tmpfld(I,J,bi,bj) + 1.
157	mlosch	1.26	ENDIF
158			ENDDO
159			kLowC(I,J,bi,bj) = 0
160			DO K= 1, Nr
161			IF (hFacC(I,J,K,bi,bj).NE.0) THEN
162			kLowC(I,J,bi,bj) = K
163	adcroft	1.21	ENDIF
164	adcroft	1.16	ENDDO
165			ENDDO
166			ENDDO
167	adcroft	1.21	C - end bi,bj loops.
168	adcroft	1.16	ENDDO
169			ENDDO
170	adcroft	1.21
171	jmc	1.32	C CALL PLOT_FIELD_XYRS( tmpfld,
172	adcroft	1.21	C & 'Model Depths K Index' , 1, myThid )
173	jmc	1.32	CALL PLOT_FIELD_XYRS(R_low,
174	adcroft	1.21	& 'Model R_low (ini_masks_etc)', 1, myThid)
175	jmc	1.32	CALL PLOT_FIELD_XYRS(Ro_surf,
176	adcroft	1.21	& 'Model Ro_surf (ini_masks_etc)', 1, myThid)
177	adcroft	1.16
178			C Calculate quantities derived from XY depth map
179	jmc	1.34	threadArea = 0. _d 0
180	adcroft	1.16	DO bj = myByLo(myThid), myByHi(myThid)
181			DO bi = myBxLo(myThid), myBxHi(myThid)
182	adcroft	1.21	DO j=1-Oly,sNy+Oly
183			DO i=1-Olx,sNx+Olx
184			C Total fluid column thickness (r_unit) :
185			c Rcolumn(i,j,bi,bj)= Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
186			tmpfld(i,j,bi,bj) = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
187			C Inverse of fluid column thickness (1/r_unit)
188			IF ( tmpfld(i,j,bi,bj) .LE. 0. ) THEN
189			recip_Rcol(i,j,bi,bj) = 0.
190	adcroft	1.16	ELSE
191	jmc	1.32	recip_Rcol(i,j,bi,bj) = 1. _d 0 / tmpfld(i,j,bi,bj)
192	adcroft	1.16	ENDIF
193			ENDDO
194			ENDDO
195	jmc	1.30	C- Compute the domain Area:
196	jmc	1.35	tileArea(bi,bj) = 0. _d 0
197	jmc	1.30	DO j=1,sNy
198			DO i=1,sNx
199	jmc	1.35	tileArea(bi,bj) = tileArea(bi,bj)
200			& + rA(i,j,bi,bj)*maskH(i,j,bi,bj)
201	jmc	1.30	ENDDO
202			ENDDO
203	jmc	1.35	threadArea = threadArea + tileArea(bi,bj)
204	adcroft	1.16	ENDDO
205			ENDDO
206	jmc	1.40	C _EXCH_XY_RS( recip_Rcol, myThid )
207	jmc	1.35	#ifdef ALLOW_AUTODIFF_TAMC
208	jmc	1.34	C_jmc: apply GLOBAL_SUM to thread-local variable (not in common block)
209	jmc	1.40	_GLOBAL_SUM_RL( threadArea, myThid )
210	jmc	1.35	#else
211			CALL GLOBAL_SUM_TILE_RL( tileArea, threadArea, myThid )
212			#endif
213	jmc	1.34	_BEGIN_MASTER( myThid )
214			globalArea = threadArea
215	jmc	1.36	C- list empty tiles:
216			msgBuf(1:1) = ' '
217			DO bj = 1,nSy
218			DO bi = 1,nSx
219			IF ( tileArea(bi,bj).EQ.0. _d 0 ) THEN
220			#ifdef ALLOW_EXCH2
221			WRITE(msgBuf,'(A,I6,A,I6,A)')
222			& 'Empty tile: #', W2_myTileList(bi), ' (bi=', bi,' )'
223			#else
224			WRITE(msgBuf,'(A,I6,I6)') 'Empty tile bi,bj=', bi, bj
225			#endif
226			CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
227			& SQUEEZE_RIGHT, myThid )
228			ENDIF
229			ENDDO
230			ENDDO
231			IF ( msgBuf(1:1).NE.' ' ) THEN
232			WRITE(msgBuf,'(A)') ' '
233			CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
234			& SQUEEZE_RIGHT, myThid )
235			ENDIF
236	jmc	1.34	_END_MASTER( myThid )
237	adcroft	1.1
238			C hFacW and hFacS (at U and V points)
239			DO bj=myByLo(myThid), myByHi(myThid)
240			DO bi=myBxLo(myThid), myBxHi(myThid)
241	cnh	1.4	DO K=1, Nr
242	adcroft	1.28	DO J=1-Oly,sNy+Oly
243	jmc	1.33	hFacW(1-OLx,J,K,bi,bj)= 0.
244	adcroft	1.28	DO I=2-Olx,sNx+Olx
245	adcroft	1.1	hFacW(I,J,K,bi,bj)=
246			& MIN(hFacC(I,J,K,bi,bj),hFacC(I-1,J,K,bi,bj))
247	adcroft	1.28	ENDDO
248			ENDDO
249	jmc	1.33	DO I=1-Olx,sNx+Olx
250			hFacS(I,1-OLy,K,bi,bj)= 0.
251			ENDDO
252	adcroft	1.28	DO J=2-Oly,sNy+oly
253			DO I=1-Olx,sNx+Olx
254	adcroft	1.1	hFacS(I,J,K,bi,bj)=
255			& MIN(hFacC(I,J,K,bi,bj),hFacC(I,J-1,K,bi,bj))
256			ENDDO
257			ENDDO
258			ENDDO
259			ENDDO
260			ENDDO
261	adcroft	1.21	CALL EXCH_UV_XYZ_RS(hFacW,hFacS,.FALSE.,myThid)
262			C The following block allows thin walls representation of non-periodic
263			C boundaries such as happen on the lat-lon grid at the N/S poles.
264			C We should really supply a flag for doing this.
265	adcroft	1.19	DO bj=myByLo(myThid), myByHi(myThid)
266			DO bi=myBxLo(myThid), myBxHi(myThid)
267			DO K=1, Nr
268			DO J=1-Oly,sNy+Oly
269	adcroft	1.21	DO I=1-Olx,sNx+Olx
270	jmc	1.34	IF (dyG(I,J,bi,bj).EQ.0.) hFacW(I,J,K,bi,bj)=0.
271			IF (dxG(I,J,bi,bj).EQ.0.) hFacS(I,J,K,bi,bj)=0.
272	adcroft	1.19	ENDDO
273			ENDDO
274			ENDDO
275			ENDDO
276			ENDDO
277	jmc	1.22
278	jmc	1.39	#ifdef ALLOW_NONHYDROSTATIC
279			C rLow & reference rSurf at Western & Southern edges (U and V points)
280			DO bj=myByLo(myThid), myByHi(myThid)
281			DO bi=myBxLo(myThid), myBxHi(myThid)
282			I = 1-OlX
283			DO J=1-Oly,sNy+Oly
284			rLowW (i,j,bi,bj) = 0.
285			rSurfW(i,j,bi,bj) = 0.
286			ENDDO
287			J = 1-Oly
288			DO I=1-Olx,sNx+Olx
289			rLowS (i,j,bi,bj) = 0.
290			rSurfS(i,j,bi,bj) = 0.
291			ENDDO
292			DO J=1-Oly,sNy+Oly
293			DO I=2-Olx,sNx+Olx
294			rSurfW(i,j,bi,bj) =
295			& MIN( Ro_surf(i-1,j,bi,bj), Ro_surf(i,j,bi,bj) )
296			rLowW(i,j,bi,bj) =
297			& MAX( R_low(i-1,j,bi,bj), R_low(i,j,bi,bj) )
298			ENDDO
299			ENDDO
300			DO J=2-Oly,sNy+Oly
301			DO I=1-Olx,sNx+Olx
302			rSurfS(i,j,bi,bj) =
303			& MIN( Ro_surf(i,j-1,bi,bj), Ro_surf(i,j,bi,bj) )
304			rLowS(i,j,bi,bj) =
305			& MAX( R_low(i,j-1,bi,bj), R_low(i,j,bi,bj) )
306			ENDDO
307			ENDDO
308			ENDDO
309			ENDDO
310			CALL EXCH_UV_XY_RS( rSurfW, rSurfS, .FALSE., myThid )
311			CALL EXCH_UV_XY_RS( rLowW, rLowS, .FALSE., myThid )
312			#endif /* ALLOW_NONHYDROSTATIC */
313
314	jmc	1.22	C- Write to disk: Total Column Thickness & hFac(C,W,S):
315			_BARRIER
316	adcroft	1.29	c _BEGIN_MASTER( myThid )
317			C This I/O is now done in write_grid.F
318	jmc	1.22	C CALL MDSWRITEFIELD( 'Depth', writeBinaryPrec, .TRUE.,
319			C & 'RS', 1, tmpfld, 1, -1, myThid )
320	adcroft	1.29	c CALL WRITE_FLD_XY_RS( 'Depth',' ',tmpfld,0,myThid)
321			c CALL WRITE_FLD_XYZ_RS( 'hFacC',' ',hFacC,0,myThid)
322			c CALL WRITE_FLD_XYZ_RS( 'hFacW',' ',hFacW,0,myThid)
323			c CALL WRITE_FLD_XYZ_RS( 'hFacS',' ',hFacS,0,myThid)
324			c _END_MASTER(myThid)
325	adcroft	1.19
326			CALL PLOT_FIELD_XYZRS( hFacC, 'hFacC' , Nr, 1, myThid )
327			CALL PLOT_FIELD_XYZRS( hFacW, 'hFacW' , Nr, 1, myThid )
328			CALL PLOT_FIELD_XYZRS( hFacS, 'hFacS' , Nr, 1, myThid )
329	adcroft	1.1
330			C Masks and reciprocals of hFac[CWS]
331			DO bj = myByLo(myThid), myByHi(myThid)
332			DO bi = myBxLo(myThid), myBxHi(myThid)
333	cnh	1.4	DO K=1,Nr
334	adcroft	1.19	DO J=1-Oly,sNy+Oly
335			DO I=1-Olx,sNx+Olx
336	jmc	1.32	IF (hFacC(I,J,K,bi,bj) .NE. 0. ) THEN
337			recip_hFacC(I,J,K,bi,bj) = 1. _d 0 / hFacC(I,J,K,bi,bj)
338	adcroft	1.21	maskC(I,J,K,bi,bj) = 1.
339	adcroft	1.1	ELSE
340	jmc	1.32	recip_hFacC(I,J,K,bi,bj) = 0.
341	adcroft	1.21	maskC(I,J,K,bi,bj) = 0.
342	adcroft	1.1	ENDIF
343	jmc	1.32	IF (hFacW(I,J,K,bi,bj) .NE. 0. ) THEN
344			recip_hFacW(I,J,K,bi,bj) = 1. _d 0 / hFacW(I,J,K,bi,bj)
345	adcroft	1.19	maskW(I,J,K,bi,bj) = 1.
346	adcroft	1.1	ELSE
347	jmc	1.32	recip_hFacW(I,J,K,bi,bj) = 0.
348	adcroft	1.19	maskW(I,J,K,bi,bj) = 0.
349	adcroft	1.1	ENDIF
350	jmc	1.32	IF (hFacS(I,J,K,bi,bj) .NE. 0. ) THEN
351			recip_hFacS(I,J,K,bi,bj) = 1. _d 0 / hFacS(I,J,K,bi,bj)
352	adcroft	1.19	maskS(I,J,K,bi,bj) = 1.
353	adcroft	1.1	ELSE
354	jmc	1.32	recip_hFacS(I,J,K,bi,bj) = 0.
355	adcroft	1.19	maskS(I,J,K,bi,bj) = 0.
356	adcroft	1.1	ENDIF
357			ENDDO
358			ENDDO
359			ENDDO
360	jmc	1.23	C- Calculate surface k index for interface W & S (U & V points)
361			DO J=1-Oly,sNy+Oly
362			DO I=1-Olx,sNx+Olx
363			ksurfW(I,J,bi,bj) = Nr+1
364			ksurfS(I,J,bi,bj) = Nr+1
365			DO k=Nr,1,-1
366			IF (hFacW(I,J,K,bi,bj).NE.0.) ksurfW(I,J,bi,bj) = k
367			IF (hFacS(I,J,K,bi,bj).NE.0.) ksurfS(I,J,bi,bj) = k
368			ENDDO
369			ENDDO
370			ENDDO
371			C - end bi,bj loops.
372	adcroft	1.1	ENDDO
373			ENDDO
374
375			C Calculate recipricols grid lengths
376			DO bj = myByLo(myThid), myByHi(myThid)
377			DO bi = myBxLo(myThid), myBxHi(myThid)
378	adcroft	1.19	DO J=1-Oly,sNy+Oly
379			DO I=1-Olx,sNx+Olx
380			IF ( dxG(I,J,bi,bj) .NE. 0. )
381	jmc	1.32	& recip_dxG(I,J,bi,bj)=1. _d 0/dxG(I,J,bi,bj)
382	adcroft	1.19	IF ( dyG(I,J,bi,bj) .NE. 0. )
383	jmc	1.32	& recip_dyG(I,J,bi,bj)=1. _d 0/dyG(I,J,bi,bj)
384	adcroft	1.19	IF ( dxC(I,J,bi,bj) .NE. 0. )
385	jmc	1.32	& recip_dxC(I,J,bi,bj)=1. _d 0/dxC(I,J,bi,bj)
386	adcroft	1.19	IF ( dyC(I,J,bi,bj) .NE. 0. )
387	jmc	1.32	& recip_dyC(I,J,bi,bj)=1. _d 0/dyC(I,J,bi,bj)
388	adcroft	1.19	IF ( dxF(I,J,bi,bj) .NE. 0. )
389	jmc	1.32	& recip_dxF(I,J,bi,bj)=1. _d 0/dxF(I,J,bi,bj)
390	adcroft	1.19	IF ( dyF(I,J,bi,bj) .NE. 0. )
391	jmc	1.32	& recip_dyF(I,J,bi,bj)=1. _d 0/dyF(I,J,bi,bj)
392	adcroft	1.19	IF ( dxV(I,J,bi,bj) .NE. 0. )
393	jmc	1.32	& recip_dxV(I,J,bi,bj)=1. _d 0/dxV(I,J,bi,bj)
394	adcroft	1.19	IF ( dyU(I,J,bi,bj) .NE. 0. )
395	jmc	1.32	& recip_dyU(I,J,bi,bj)=1. _d 0/dyU(I,J,bi,bj)
396	adcroft	1.19	IF ( rA(I,J,bi,bj) .NE. 0. )
397	jmc	1.32	& recip_rA(I,J,bi,bj)=1. _d 0/rA(I,J,bi,bj)
398	adcroft	1.19	IF ( rAs(I,J,bi,bj) .NE. 0. )
399	jmc	1.32	& recip_rAs(I,J,bi,bj)=1. _d 0/rAs(I,J,bi,bj)
400	adcroft	1.19	IF ( rAw(I,J,bi,bj) .NE. 0. )
401	jmc	1.32	& recip_rAw(I,J,bi,bj)=1. _d 0/rAw(I,J,bi,bj)
402	adcroft	1.19	IF ( rAz(I,J,bi,bj) .NE. 0. )
403	jmc	1.32	& recip_rAz(I,J,bi,bj)=1. _d 0/rAz(I,J,bi,bj)
404	adcroft	1.1	ENDDO
405			ENDDO
406			ENDDO
407			ENDDO
408
409	jmc	1.32	c #ifdef ALLOW_NONHYDROSTATIC
410			C-- Calculate "recip_hFacU" = reciprocal hfac distance/volume for W cells
411	jmc	1.33	C NOTE: not used ; computed locally in CALC_GW
412	jmc	1.32	c #endif
413
414	adcroft	1.1	RETURN
415			END