model/src/ini_masks_etc.F

C $Header: /u/gcmpack/MITgcm/model/src/ini_masks_etc.F,v 1.37 2008/09/10 09:19:22 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_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_R4(   recip_Rcol, myThid )
#ifdef ALLOW_AUTODIFF_TAMC
C_jmc: apply GLOBAL_SUM to thread-local variable (not in common block)
      _GLOBAL_SUM_R8( 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

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	mlosch	1.38	C $Header: /u/gcmpack/MITgcm/model/src/ini_masks_etc.F,v 1.37 2008/09/10 09:19:22 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	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			CALL SHELFICE_UPDATE_MASKS(
133			I rF, recip_drF,
134			U hFacC,
135			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	adcroft	1.21	C _EXCH_XY_R4( 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			_GLOBAL_SUM_R8( 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			C- Write to disk: Total Column Thickness & hFac(C,W,S):
279			_BARRIER
280	adcroft	1.29	c _BEGIN_MASTER( myThid )
281			C This I/O is now done in write_grid.F
282	jmc	1.22	C CALL MDSWRITEFIELD( 'Depth', writeBinaryPrec, .TRUE.,
283			C & 'RS', 1, tmpfld, 1, -1, myThid )
284	adcroft	1.29	c CALL WRITE_FLD_XY_RS( 'Depth',' ',tmpfld,0,myThid)
285			c CALL WRITE_FLD_XYZ_RS( 'hFacC',' ',hFacC,0,myThid)
286			c CALL WRITE_FLD_XYZ_RS( 'hFacW',' ',hFacW,0,myThid)
287			c CALL WRITE_FLD_XYZ_RS( 'hFacS',' ',hFacS,0,myThid)
288			c _END_MASTER(myThid)
289	adcroft	1.19
290			CALL PLOT_FIELD_XYZRS( hFacC, 'hFacC' , Nr, 1, myThid )
291			CALL PLOT_FIELD_XYZRS( hFacW, 'hFacW' , Nr, 1, myThid )
292			CALL PLOT_FIELD_XYZRS( hFacS, 'hFacS' , Nr, 1, myThid )
293	adcroft	1.1
294			C Masks and reciprocals of hFac[CWS]
295			DO bj = myByLo(myThid), myByHi(myThid)
296			DO bi = myBxLo(myThid), myBxHi(myThid)
297	cnh	1.4	DO K=1,Nr
298	adcroft	1.19	DO J=1-Oly,sNy+Oly
299			DO I=1-Olx,sNx+Olx
300	jmc	1.32	IF (hFacC(I,J,K,bi,bj) .NE. 0. ) THEN
301			recip_hFacC(I,J,K,bi,bj) = 1. _d 0 / hFacC(I,J,K,bi,bj)
302	adcroft	1.21	maskC(I,J,K,bi,bj) = 1.
303	adcroft	1.1	ELSE
304	jmc	1.32	recip_hFacC(I,J,K,bi,bj) = 0.
305	adcroft	1.21	maskC(I,J,K,bi,bj) = 0.
306	adcroft	1.1	ENDIF
307	jmc	1.32	IF (hFacW(I,J,K,bi,bj) .NE. 0. ) THEN
308			recip_hFacW(I,J,K,bi,bj) = 1. _d 0 / hFacW(I,J,K,bi,bj)
309	adcroft	1.19	maskW(I,J,K,bi,bj) = 1.
310	adcroft	1.1	ELSE
311	jmc	1.32	recip_hFacW(I,J,K,bi,bj) = 0.
312	adcroft	1.19	maskW(I,J,K,bi,bj) = 0.
313	adcroft	1.1	ENDIF
314	jmc	1.32	IF (hFacS(I,J,K,bi,bj) .NE. 0. ) THEN
315			recip_hFacS(I,J,K,bi,bj) = 1. _d 0 / hFacS(I,J,K,bi,bj)
316	adcroft	1.19	maskS(I,J,K,bi,bj) = 1.
317	adcroft	1.1	ELSE
318	jmc	1.32	recip_hFacS(I,J,K,bi,bj) = 0.
319	adcroft	1.19	maskS(I,J,K,bi,bj) = 0.
320	adcroft	1.1	ENDIF
321			ENDDO
322			ENDDO
323			ENDDO
324	jmc	1.23	C- Calculate surface k index for interface W & S (U & V points)
325			DO J=1-Oly,sNy+Oly
326			DO I=1-Olx,sNx+Olx
327			ksurfW(I,J,bi,bj) = Nr+1
328			ksurfS(I,J,bi,bj) = Nr+1
329			DO k=Nr,1,-1
330			IF (hFacW(I,J,K,bi,bj).NE.0.) ksurfW(I,J,bi,bj) = k
331			IF (hFacS(I,J,K,bi,bj).NE.0.) ksurfS(I,J,bi,bj) = k
332			ENDDO
333			ENDDO
334			ENDDO
335			C - end bi,bj loops.
336	adcroft	1.1	ENDDO
337			ENDDO
338
339			C Calculate recipricols grid lengths
340			DO bj = myByLo(myThid), myByHi(myThid)
341			DO bi = myBxLo(myThid), myBxHi(myThid)
342	adcroft	1.19	DO J=1-Oly,sNy+Oly
343			DO I=1-Olx,sNx+Olx
344			IF ( dxG(I,J,bi,bj) .NE. 0. )
345	jmc	1.32	& recip_dxG(I,J,bi,bj)=1. _d 0/dxG(I,J,bi,bj)
346	adcroft	1.19	IF ( dyG(I,J,bi,bj) .NE. 0. )
347	jmc	1.32	& recip_dyG(I,J,bi,bj)=1. _d 0/dyG(I,J,bi,bj)
348	adcroft	1.19	IF ( dxC(I,J,bi,bj) .NE. 0. )
349	jmc	1.32	& recip_dxC(I,J,bi,bj)=1. _d 0/dxC(I,J,bi,bj)
350	adcroft	1.19	IF ( dyC(I,J,bi,bj) .NE. 0. )
351	jmc	1.32	& recip_dyC(I,J,bi,bj)=1. _d 0/dyC(I,J,bi,bj)
352	adcroft	1.19	IF ( dxF(I,J,bi,bj) .NE. 0. )
353	jmc	1.32	& recip_dxF(I,J,bi,bj)=1. _d 0/dxF(I,J,bi,bj)
354	adcroft	1.19	IF ( dyF(I,J,bi,bj) .NE. 0. )
355	jmc	1.32	& recip_dyF(I,J,bi,bj)=1. _d 0/dyF(I,J,bi,bj)
356	adcroft	1.19	IF ( dxV(I,J,bi,bj) .NE. 0. )
357	jmc	1.32	& recip_dxV(I,J,bi,bj)=1. _d 0/dxV(I,J,bi,bj)
358	adcroft	1.19	IF ( dyU(I,J,bi,bj) .NE. 0. )
359	jmc	1.32	& recip_dyU(I,J,bi,bj)=1. _d 0/dyU(I,J,bi,bj)
360	adcroft	1.19	IF ( rA(I,J,bi,bj) .NE. 0. )
361	jmc	1.32	& recip_rA(I,J,bi,bj)=1. _d 0/rA(I,J,bi,bj)
362	adcroft	1.19	IF ( rAs(I,J,bi,bj) .NE. 0. )
363	jmc	1.32	& recip_rAs(I,J,bi,bj)=1. _d 0/rAs(I,J,bi,bj)
364	adcroft	1.19	IF ( rAw(I,J,bi,bj) .NE. 0. )
365	jmc	1.32	& recip_rAw(I,J,bi,bj)=1. _d 0/rAw(I,J,bi,bj)
366	adcroft	1.19	IF ( rAz(I,J,bi,bj) .NE. 0. )
367	jmc	1.32	& recip_rAz(I,J,bi,bj)=1. _d 0/rAz(I,J,bi,bj)
368	adcroft	1.1	ENDDO
369			ENDDO
370			ENDDO
371			ENDDO
372
373	jmc	1.32	c #ifdef ALLOW_NONHYDROSTATIC
374			C-- Calculate "recip_hFacU" = reciprocal hfac distance/volume for W cells
375	jmc	1.33	C NOTE: not used ; computed locally in CALC_GW
376	jmc	1.32	c #endif
377
378	adcroft	1.1	RETURN
379			END