model/src/ini_masks_etc.F

C $Header: /u/gcmpack/models/MITgcmUV/model/src/ini_masks_etc.F,v 1.16 1999/05/05 18:32:34 adcroft Exp $

#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"

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

C     Calculate lopping factor hFacC
      DO bj=myByLo(myThid), myByHi(myThid)
       DO bi=myBxLo(myThid), myBxHi(myThid)
        DO K=1, Nr
         DO J=1,sNy
          DO I=1,sNx
C          Round depths within a small fraction of layer depth to that
C          layer depth.
           IF ( ABS(H(I,J,bi,bj)-rF(K)) .LT. 
     &          1. _d -6*ABS(rF(K)) .AND.
     &          ABS(H(I,J,bi,bj)-rF(K)) .LT. 
     &          1. _d -6*ABS(H(I,J,bi,bj)) )THEN
            H(I,J,bi,bj) = rF(K)
           ENDIF
           IF     ( H(I,J,bi,bj)*rkFac .GE. rF(K)*rkFac ) THEN
C           Top of cell is below base of domain
            hFacC(I,J,K,bi,bj) = 0.
           ELSEIF ( H(I,J,bi,bj)*rkFac .LE. rF(K+1)*rkFac ) THEN
C           Base of domain is below bottom of this cell
            hFacC(I,J,K,bi,bj) = 1.
           ELSE
C           Base of domain is in this cell
C           Set hFac to the fraction of the cell that is open.
            hFacC(I,J,K,bi,bj) =
     &       (rF(K)*rkFac-H(I,J,bi,bj)*rkFac)*recip_drF(K)
           ENDIF
C Impose minimum fraction
           IF (hFacC(I,J,K,bi,bj).LT.hFacMin) THEN
            IF (hFacC(I,J,K,bi,bj).LT.hFacMin*0.5) THEN
             hFacC(I,J,K,bi,bj)=0.
            ELSE
             hFacC(I,J,K,bi,bj)=hFacMin
            ENDIF
           ENDIF
C Impose minimum size (dimensional)
           IF (drF(k)*hFacC(I,J,K,bi,bj).LT.hFacMinDr) THEN
            IF (drF(k)*hFacC(I,J,K,bi,bj).LT.hFacMinDr*0.5) THEN
             hFacC(I,J,K,bi,bj)=0.
            ELSE
             hFacC(I,J,K,bi,bj)=hFacMinDr*recip_drF(k)
            ENDIF
           ENDIF
           depthInK(i,j,bi,bj) = depthInK(i,j,bi,bj) + 1.
Crg  &                          +hFacC(i,j,k,bi,bj)
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      _EXCH_XYZ_R4(hFacC , myThid )
      _EXCH_XY_R4( depthInK, myThid )

      CALL PLOT_FIELD_XYRS( depthInK, 
     & 'Model Depths K Index' , 1, myThid )

C Re-calculate depth of ocean, taking into account hFacC
      DO bj=myByLo(myThid), myByHi(myThid)
       DO bi=myBxLo(myThid), myBxHi(myThid)
        DO J=1,sNy
         DO I=1,sNx
          H(I,J,bi,bj)=0.
          DO K=1,Nr
           H(I,J,bi,bj)=H(I,J,bi,bj)-
     &       rkFac*drF(k)*hFacC(I,J,K,bi,bj)
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      _EXCH_XY_R4(H , myThid )
      CALL WRITE_FLD_XY_RS( 'Depth',' ',H,0,myThid)
C     CALL MDSWRITEFIELD( 'Depth', writeBinaryPrec, .TRUE.,
C    &                    'RS', 1, H, 1, -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,sNy
         DO I=1,sNx
C         Inverse of depth
          IF ( h(i,j,bi,bj) .EQ. 0. _d 0 ) THEN
           recip_H(i,j,bi,bj) = 0. _d 0
          ELSE
           recip_H(i,j,bi,bj) = 1. _d 0 /  abs( H(i,j,bi,bj) )
          ENDIF
          depthInK(i,j,bi,bj) = 0.
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      _EXCH_XY_R4(   recip_H, 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
      _EXCH_XYZ_R4(hFacW , myThid )
      _EXCH_XYZ_R4(hFacS , 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,sNy
          DO I=1,sNx
           IF (HFacC(I,J,K,bi,bj) .NE. 0. _d 0 ) THEN
            recip_HFacC(I,J,K,bi,bj) = 1. _d 0 / HFacC(I,J,K,bi,bj)
           ELSE
            recip_HFacC(I,J,K,bi,bj) = 0. _d 0
           ENDIF
           IF (HFacW(I,J,K,bi,bj) .NE. 0. _d 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. _d 0
           ELSE
            recip_HFacW(I,J,K,bi,bj) = 0. _d 0
            maskW(I,J,K,bi,bj) = 0.0 _d 0
           ENDIF
           IF (HFacS(I,J,K,bi,bj) .NE. 0. _d 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. _d 0
           ELSE
            recip_HFacS(I,J,K,bi,bj) = 0. _d 0
            maskS(I,J,K,bi,bj) = 0. _d 0
           ENDIF
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      _EXCH_XYZ_R4(recip_HFacC    , myThid )
      _EXCH_XYZ_R4(recip_HFacW    , myThid )
      _EXCH_XYZ_R4(recip_HFacS    , myThid )
      _EXCH_XYZ_R4(maskW    , myThid )
      _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,sNy
         DO I=1,sNx
          recip_dxG(I,J,bi,bj)=1.d0/dxG(I,J,bi,bj)
          recip_dyG(I,J,bi,bj)=1.d0/dyG(I,J,bi,bj)
          recip_dxC(I,J,bi,bj)=1.d0/dxC(I,J,bi,bj)
          recip_dyC(I,J,bi,bj)=1.d0/dyC(I,J,bi,bj)
          recip_dxF(I,J,bi,bj)=1.d0/dxF(I,J,bi,bj)
          recip_dyF(I,J,bi,bj)=1.d0/dyF(I,J,bi,bj)
          recip_dxV(I,J,bi,bj)=1.d0/dxV(I,J,bi,bj)
          recip_dyU(I,J,bi,bj)=1.d0/dyU(I,J,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      _EXCH_XY_R4(recip_dxG, myThid )
      _EXCH_XY_R4(recip_dyG, myThid )
      _EXCH_XY_R4(recip_dxC, myThid )
      _EXCH_XY_R4(recip_dyC, myThid )
      _EXCH_XY_R4(recip_dxF, myThid )
      _EXCH_XY_R4(recip_dyF, myThid )
      _EXCH_XY_R4(recip_dxV, myThid )
      _EXCH_XY_R4(recip_dyU, 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,sNy
          DO I=1,sNx
           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
      _EXCH_XY_R4(recip_hFacU, myThid )
#endif
C
      RETURN
      END
1	adcroft	1.17	C $Header: /u/gcmpack/models/MITgcmUV/model/src/ini_masks_etc.F,v 1.16 1999/05/05 18:32:34 adcroft Exp $
2	adcroft	1.1
3	cnh	1.11	#include "CPP_OPTIONS.h"
4	adcroft	1.1
5			CStartOfInterface
6			SUBROUTINE INI_MASKS_ETC( myThid )
7			C /==========================================================\
8			C \| SUBROUTINE INI_MASKS_ETC \|
9			C \| o Initialise masks and topography factors \|
10			C \|==========================================================\|
11			C \| These arrays are used throughout the code and describe \|
12			C \| the topography of the domain through masks (0s and 1s) \|
13			C \| and fractional height factors (0<hFac<1). The latter \|
14			C \| distinguish between the lopped-cell and full-step \|
15			C \| topographic representations. \|
16			C \==========================================================/
17	adcroft	1.13	IMPLICIT NONE
18	adcroft	1.1
19			C === Global variables ===
20			#include "SIZE.h"
21			#include "EEPARAMS.h"
22			#include "PARAMS.h"
23			#include "GRID.h"
24
25			C == Routine arguments ==
26	cnh	1.6	C myThid - Number of this instance of INI_MASKS_ETC
27	adcroft	1.1	INTEGER myThid
28			CEndOfInterface
29
30			C == Local variables ==
31			C bi,bj - Loop counters
32			C I,J,K
33			INTEGER bi, bj
34			INTEGER I, J, K
35	adcroft	1.15	#ifdef ALLOW_NONHYDROSTATIC
36			INTEGER Km1
37			_RL hFacUpper,hFacLower
38			#endif
39	adcroft	1.1
40	adcroft	1.2	C Calculate lopping factor hFacC
41			DO bj=myByLo(myThid), myByHi(myThid)
42			DO bi=myBxLo(myThid), myBxHi(myThid)
43	cnh	1.4	DO K=1, Nr
44	adcroft	1.2	DO J=1,sNy
45			DO I=1,sNx
46	cnh	1.7	C Round depths within a small fraction of layer depth to that
47			C layer depth.
48	cnh	1.9	IF ( ABS(H(I,J,bi,bj)-rF(K)) .LT.
49			& 1. _d -6*ABS(rF(K)) .AND.
50			& ABS(H(I,J,bi,bj)-rF(K)) .LT.
51			& 1. _d -6*ABS(H(I,J,bi,bj)) )THEN
52	cnh	1.7	H(I,J,bi,bj) = rF(K)
53			ENDIF
54	cnh	1.6	IF ( H(I,J,bi,bj)rkFac .GE. rF(K)rkFac ) THEN
55	adcroft	1.3	C Top of cell is below base of domain
56	adcroft	1.2	hFacC(I,J,K,bi,bj) = 0.
57	cnh	1.6	ELSEIF ( H(I,J,bi,bj)rkFac .LE. rF(K+1)rkFac ) THEN
58	adcroft	1.3	C Base of domain is below bottom of this cell
59	adcroft	1.2	hFacC(I,J,K,bi,bj) = 1.
60			ELSE
61			C Base of domain is in this cell
62	adcroft	1.3	C Set hFac to the fraction of the cell that is open.
63	adcroft	1.15	hFacC(I,J,K,bi,bj) =
64			& (rF(K)rkFac-H(I,J,bi,bj)rkFac)*recip_drF(K)
65	adcroft	1.3	ENDIF
66			C Impose minimum fraction
67			IF (hFacC(I,J,K,bi,bj).LT.hFacMin) THEN
68			IF (hFacC(I,J,K,bi,bj).LT.hFacMin*0.5) THEN
69			hFacC(I,J,K,bi,bj)=0.
70			ELSE
71			hFacC(I,J,K,bi,bj)=hFacMin
72			ENDIF
73			ENDIF
74			C Impose minimum size (dimensional)
75	cnh	1.5	IF (drF(k)*hFacC(I,J,K,bi,bj).LT.hFacMinDr) THEN
76			IF (drF(k)hFacC(I,J,K,bi,bj).LT.hFacMinDr0.5) THEN
77	adcroft	1.3	hFacC(I,J,K,bi,bj)=0.
78			ELSE
79	cnh	1.5	hFacC(I,J,K,bi,bj)=hFacMinDr*recip_drF(k)
80	adcroft	1.3	ENDIF
81	adcroft	1.2	ENDIF
82	adcroft	1.17	depthInK(i,j,bi,bj) = depthInK(i,j,bi,bj) + 1.
83			Crg & +hFacC(i,j,k,bi,bj)
84	adcroft	1.2	ENDDO
85			ENDDO
86			ENDDO
87			ENDDO
88			ENDDO
89			_EXCH_XYZ_R4(hFacC , myThid )
90	cnh	1.7	_EXCH_XY_R4( depthInK, myThid )
91
92	cnh	1.9	CALL PLOT_FIELD_XYRS( depthInK,
93			& 'Model Depths K Index' , 1, myThid )
94	adcroft	1.16
95			C Re-calculate depth of ocean, taking into account hFacC
96			DO bj=myByLo(myThid), myByHi(myThid)
97			DO bi=myBxLo(myThid), myBxHi(myThid)
98			DO J=1,sNy
99			DO I=1,sNx
100			H(I,J,bi,bj)=0.
101			DO K=1,Nr
102			H(I,J,bi,bj)=H(I,J,bi,bj)-
103			& rkFacdrF(k)hFacC(I,J,K,bi,bj)
104			ENDDO
105			ENDDO
106			ENDDO
107			ENDDO
108			ENDDO
109			_EXCH_XY_R4(H , myThid )
110			CALL WRITE_FLD_XY_RS( 'Depth',' ',H,0,myThid)
111			C CALL MDSWRITEFIELD( 'Depth', writeBinaryPrec, .TRUE.,
112			C & 'RS', 1, H, 1, -1, myThid )
113
114			C Calculate quantities derived from XY depth map
115			DO bj = myByLo(myThid), myByHi(myThid)
116			DO bi = myBxLo(myThid), myBxHi(myThid)
117			DO J=1,sNy
118			DO I=1,sNx
119			C Inverse of depth
120			IF ( h(i,j,bi,bj) .EQ. 0. _d 0 ) THEN
121			recip_H(i,j,bi,bj) = 0. _d 0
122			ELSE
123			recip_H(i,j,bi,bj) = 1. _d 0 / abs( H(i,j,bi,bj) )
124			ENDIF
125			depthInK(i,j,bi,bj) = 0.
126			ENDDO
127			ENDDO
128			ENDDO
129			ENDDO
130			_EXCH_XY_R4( recip_H, myThid )
131	adcroft	1.1
132			C hFacW and hFacS (at U and V points)
133			DO bj=myByLo(myThid), myByHi(myThid)
134			DO bi=myBxLo(myThid), myBxHi(myThid)
135	cnh	1.4	DO K=1, Nr
136	adcroft	1.1	DO J=1,sNy
137			DO I=1,sNx
138			hFacW(I,J,K,bi,bj)=
139			& MIN(hFacC(I,J,K,bi,bj),hFacC(I-1,J,K,bi,bj))
140			hFacS(I,J,K,bi,bj)=
141			& MIN(hFacC(I,J,K,bi,bj),hFacC(I,J-1,K,bi,bj))
142			ENDDO
143			ENDDO
144			ENDDO
145			ENDDO
146			ENDDO
147			_EXCH_XYZ_R4(hFacW , myThid )
148			_EXCH_XYZ_R4(hFacS , myThid )
149
150			C Masks and reciprocals of hFac[CWS]
151			DO bj = myByLo(myThid), myByHi(myThid)
152			DO bi = myBxLo(myThid), myBxHi(myThid)
153	cnh	1.4	DO K=1,Nr
154	adcroft	1.1	DO J=1,sNy
155			DO I=1,sNx
156	cnh	1.10	IF (HFacC(I,J,K,bi,bj) .NE. 0. _d 0 ) THEN
157			recip_HFacC(I,J,K,bi,bj) = 1. _d 0 / HFacC(I,J,K,bi,bj)
158	adcroft	1.1	ELSE
159	cnh	1.10	recip_HFacC(I,J,K,bi,bj) = 0. _d 0
160	adcroft	1.1	ENDIF
161	cnh	1.10	IF (HFacW(I,J,K,bi,bj) .NE. 0. _d 0 ) THEN
162			recip_HFacW(I,J,K,bi,bj) = 1. _d 0 / HFacW(I,J,K,bi,bj)
163			maskW(I,J,K,bi,bj) = 1. _d 0
164	adcroft	1.1	ELSE
165	cnh	1.10	recip_HFacW(I,J,K,bi,bj) = 0. _d 0
166			maskW(I,J,K,bi,bj) = 0.0 _d 0
167	adcroft	1.1	ENDIF
168	cnh	1.10	IF (HFacS(I,J,K,bi,bj) .NE. 0. _d 0 ) THEN
169			recip_HFacS(I,J,K,bi,bj) = 1. _d 0 / HFacS(I,J,K,bi,bj)
170			maskS(I,J,K,bi,bj) = 1. _d 0
171	adcroft	1.1	ELSE
172	cnh	1.10	recip_HFacS(I,J,K,bi,bj) = 0. _d 0
173			maskS(I,J,K,bi,bj) = 0. _d 0
174	adcroft	1.1	ENDIF
175			ENDDO
176			ENDDO
177			ENDDO
178			ENDDO
179			ENDDO
180	cnh	1.4	_EXCH_XYZ_R4(recip_HFacC , myThid )
181			_EXCH_XYZ_R4(recip_HFacW , myThid )
182			_EXCH_XYZ_R4(recip_HFacS , myThid )
183	adcroft	1.1	_EXCH_XYZ_R4(maskW , myThid )
184			_EXCH_XYZ_R4(maskS , myThid )
185
186			C Calculate recipricols grid lengths
187			DO bj = myByLo(myThid), myByHi(myThid)
188			DO bi = myBxLo(myThid), myBxHi(myThid)
189			DO J=1,sNy
190			DO I=1,sNx
191	cnh	1.4	recip_dxG(I,J,bi,bj)=1.d0/dxG(I,J,bi,bj)
192			recip_dyG(I,J,bi,bj)=1.d0/dyG(I,J,bi,bj)
193			recip_dxC(I,J,bi,bj)=1.d0/dxC(I,J,bi,bj)
194			recip_dyC(I,J,bi,bj)=1.d0/dyC(I,J,bi,bj)
195			recip_dxF(I,J,bi,bj)=1.d0/dxF(I,J,bi,bj)
196			recip_dyF(I,J,bi,bj)=1.d0/dyF(I,J,bi,bj)
197			recip_dxV(I,J,bi,bj)=1.d0/dxV(I,J,bi,bj)
198			recip_dyU(I,J,bi,bj)=1.d0/dyU(I,J,bi,bj)
199	adcroft	1.1	ENDDO
200			ENDDO
201			ENDDO
202			ENDDO
203	cnh	1.4	_EXCH_XY_R4(recip_dxG, myThid )
204			_EXCH_XY_R4(recip_dyG, myThid )
205			_EXCH_XY_R4(recip_dxC, myThid )
206			_EXCH_XY_R4(recip_dyC, myThid )
207			_EXCH_XY_R4(recip_dxF, myThid )
208			_EXCH_XY_R4(recip_dyF, myThid )
209			_EXCH_XY_R4(recip_dxV, myThid )
210			_EXCH_XY_R4(recip_dyU, myThid )
211	adcroft	1.1
212	adcroft	1.15	#ifdef ALLOW_NONHYDROSTATIC
213			C-- Calculate the reciprocal hfac distance/volume for W cells
214			DO bj = myByLo(myThid), myByHi(myThid)
215			DO bi = myBxLo(myThid), myBxHi(myThid)
216			DO K=1,Nr
217			Km1=max(K-1,1)
218			hFacUpper=drF(Km1)/(drF(Km1)+drF(K))
219			IF (Km1.EQ.K) hFacUpper=0.
220			hFacLower=drF(K)/(drF(Km1)+drF(K))
221			DO J=1,sNy
222			DO I=1,sNx
223			IF (hFacC(I,J,K,bi,bj).NE.0.) THEN
224			IF (hFacC(I,J,K,bi,bj).LE.0.5) THEN
225			recip_hFacU(I,J,K,bi,bj)=
226			& hFacUpper+hFacLower*hFacC(I,J,K,bi,bj)
227			ELSE
228			recip_hFacU(I,J,K,bi,bj)=1.
229			ENDIF
230			ELSE
231			recip_hFacU(I,J,K,bi,bj)=0.
232			ENDIF
233			IF (recip_hFacU(I,J,K,bi,bj).NE.0.)
234			& recip_hFacU(I,J,K,bi,bj)=1./recip_hFacU(I,J,K,bi,bj)
235			ENDDO
236			ENDDO
237			ENDDO
238			ENDDO
239			ENDDO
240			_EXCH_XY_R4(recip_hFacU, myThid )
241			#endif
242	adcroft	1.1	C
243			RETURN
244			END