model/src/ini_masks_etc.F

C $Header: /u/gcmpack/models/MITgcmUV/model/src/ini_masks_etc.F,v 1.8 1998/09/29 18:50:57 cnh Exp $

#include "CPP_EEOPTIONS.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     \==========================================================/

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

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 )
      IF ( myThid .EQ. 1 ) WRITE(0,*) 'AAAA'

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.
C           hFacC(I,J,K,bi,bj) = (rF(K)*rkFac-H(I,J,bi,bj)*rkFac)*recip_drF(K)
CcnhDebugStarts
C Impose full-step
            hFacC(I,J,K,bi,bj) = 1.
CCnhDebugEnds
           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) 
     &                          +hFacC(i,j,k,bi,bj)
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      IF ( myThid .EQ. 1 ) WRITE(0,*) 'BBBB'
      _EXCH_XYZ_R4(hFacC , myThid )
      IF ( myThid .EQ. 1 ) WRITE(0,*) 'CCCC'
      _EXCH_XY_R4( depthInK, myThid )

      IF ( myThid .EQ. 1 ) WRITE(0,*) 'DDDD'
      CALL PLOT_FIELD_XYRS( depthInK, 
     & 'Model Depths K Index' , 1, myThid )

C     hFacW and hFacS (at U and V points)
      DO bj=myByLo(myThid), myByHi(myThid)
       DO bi=myBxLo(myThid), myBxHi(myThid)
        DO K=1, Nr
         DO J=1,sNy
          DO I=1,sNx
           hFacW(I,J,K,bi,bj)=
     &       MIN(hFacC(I,J,K,bi,bj),hFacC(I-1,J,K,bi,bj))
           hFacS(I,J,K,bi,bj)=
     &       MIN(hFacC(I,J,K,bi,bj),hFacC(I,J-1,K,bi,bj))
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      _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. D0 ) THEN
            recip_HFacC(I,J,K,bi,bj) = 1. D0 / HFacC(I,J,K,bi,bj)
           ELSE
            recip_HFacC(I,J,K,bi,bj) = 0. D0
           ENDIF
           IF (HFacW(I,J,K,bi,bj) .NE. 0. D0 ) THEN
            recip_HFacW(I,J,K,bi,bj) = 1. D0 / HFacW(I,J,K,bi,bj)
            maskW(I,J,K,bi,bj) = 1. D0
           ELSE
            recip_HFacW(I,J,K,bi,bj) = 0. D0
            maskW(I,J,K,bi,bj) = 0.0 D0
           ENDIF
           IF (HFacS(I,J,K,bi,bj) .NE. 0. D0 ) THEN
            recip_HFacS(I,J,K,bi,bj) = 1. D0 / HFacS(I,J,K,bi,bj)
            maskS(I,J,K,bi,bj) = 1. D0
           ELSE
            recip_HFacS(I,J,K,bi,bj) = 0. D0
            maskS(I,J,K,bi,bj) = 0. D0
           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 )

C
      RETURN
      END
1	C $Header: /u/gcmpack/models/MITgcmUV/model/src/ini_masks_etc.F,v 1.8 1998/09/29 18:50:57 cnh Exp $
2
3	#include "CPP_EEOPTIONS.h"
4
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
18	C === Global variables ===
19	#include "SIZE.h"
20	#include "EEPARAMS.h"
21	#include "PARAMS.h"
22	#include "GRID.h"
23
24	C == Routine arguments ==
25	C myThid - Number of this instance of INI_MASKS_ETC
26	INTEGER myThid
27	CEndOfInterface
28
29	C == Local variables ==
30	C bi,bj - Loop counters
31	C I,J,K
32	INTEGER bi, bj
33	INTEGER I, J, K
34
35	C Calculate quantities derived from XY depth map
36	DO bj = myByLo(myThid), myByHi(myThid)
37	DO bi = myBxLo(myThid), myBxHi(myThid)
38	DO J=1,sNy
39	DO I=1,sNx
40	C Inverse of depth
41	IF ( h(i,j,bi,bj) .EQ. 0. _d 0 ) THEN
42	recip_H(i,j,bi,bj) = 0. _d 0
43	ELSE
44	recip_H(i,j,bi,bj) = 1. _d 0 / abs( H(i,j,bi,bj) )
45	ENDIF
46	depthInK(i,j,bi,bj) = 0.
47	ENDDO
48	ENDDO
49	ENDDO
50	ENDDO
51	_EXCH_XY_R4( recip_H, myThid )
52	IF ( myThid .EQ. 1 ) WRITE(0,*) 'AAAA'
53
54	C Calculate lopping factor hFacC
55	DO bj=myByLo(myThid), myByHi(myThid)
56	DO bi=myBxLo(myThid), myBxHi(myThid)
57	DO K=1, Nr
58	DO J=1,sNy
59	DO I=1,sNx
60	C Round depths within a small fraction of layer depth to that
61	C layer depth.
62	IF ( ABS(H(I,J,bi,bj)-rF(K)) .LT.
63	& 1. _d -6*ABS(rF(K)) .AND.
64	& ABS(H(I,J,bi,bj)-rF(K)) .LT.
65	& 1. _d -6*ABS(H(I,J,bi,bj)) )THEN
66	H(I,J,bi,bj) = rF(K)
67	ENDIF
68	IF ( H(I,J,bi,bj)rkFac .GE. rF(K)rkFac ) THEN
69	C Top of cell is below base of domain
70	hFacC(I,J,K,bi,bj) = 0.
71	ELSEIF ( H(I,J,bi,bj)rkFac .LE. rF(K+1)rkFac ) THEN
72	C Base of domain is below bottom of this cell
73	hFacC(I,J,K,bi,bj) = 1.
74	ELSE
75	C Base of domain is in this cell
76	C Set hFac to the fraction of the cell that is open.
77	C hFacC(I,J,K,bi,bj) = (rF(K)rkFac-H(I,J,bi,bj)rkFac)*recip_drF(K)
78	CcnhDebugStarts
79	C Impose full-step
80	hFacC(I,J,K,bi,bj) = 1.
81	CCnhDebugEnds
82	ENDIF
83	C Impose minimum fraction
84	IF (hFacC(I,J,K,bi,bj).LT.hFacMin) THEN
85	IF (hFacC(I,J,K,bi,bj).LT.hFacMin*0.5) THEN
86	hFacC(I,J,K,bi,bj)=0.
87	ELSE
88	hFacC(I,J,K,bi,bj)=hFacMin
89	ENDIF
90	ENDIF
91	C Impose minimum size (dimensional)
92	IF (drF(k)*hFacC(I,J,K,bi,bj).LT.hFacMinDr) THEN
93	IF (drF(k)hFacC(I,J,K,bi,bj).LT.hFacMinDr0.5) THEN
94	hFacC(I,J,K,bi,bj)=0.
95	ELSE
96	hFacC(I,J,K,bi,bj)=hFacMinDr*recip_drF(k)
97	ENDIF
98	ENDIF
99	depthInK(i,j,bi,bj) = depthInK(i,j,bi,bj)
100	& +hFacC(i,j,k,bi,bj)
101	ENDDO
102	ENDDO
103	ENDDO
104	ENDDO
105	ENDDO
106	IF ( myThid .EQ. 1 ) WRITE(0,*) 'BBBB'
107	_EXCH_XYZ_R4(hFacC , myThid )
108	IF ( myThid .EQ. 1 ) WRITE(0,*) 'CCCC'
109	_EXCH_XY_R4( depthInK, myThid )
110
111	IF ( myThid .EQ. 1 ) WRITE(0,*) 'DDDD'
112	CALL PLOT_FIELD_XYRS( depthInK,
113	& 'Model Depths K Index' , 1, myThid )
114
115	C hFacW and hFacS (at U and V points)
116	DO bj=myByLo(myThid), myByHi(myThid)
117	DO bi=myBxLo(myThid), myBxHi(myThid)
118	DO K=1, Nr
119	DO J=1,sNy
120	DO I=1,sNx
121	hFacW(I,J,K,bi,bj)=
122	& MIN(hFacC(I,J,K,bi,bj),hFacC(I-1,J,K,bi,bj))
123	hFacS(I,J,K,bi,bj)=
124	& MIN(hFacC(I,J,K,bi,bj),hFacC(I,J-1,K,bi,bj))
125	ENDDO
126	ENDDO
127	ENDDO
128	ENDDO
129	ENDDO
130	_EXCH_XYZ_R4(hFacW , myThid )
131	_EXCH_XYZ_R4(hFacS , myThid )
132
133	C Masks and reciprocals of hFac[CWS]
134	DO bj = myByLo(myThid), myByHi(myThid)
135	DO bi = myBxLo(myThid), myBxHi(myThid)
136	DO K=1,Nr
137	DO J=1,sNy
138	DO I=1,sNx
139	IF (HFacC(I,J,K,bi,bj) .NE. 0. D0 ) THEN
140	recip_HFacC(I,J,K,bi,bj) = 1. D0 / HFacC(I,J,K,bi,bj)
141	ELSE
142	recip_HFacC(I,J,K,bi,bj) = 0. D0
143	ENDIF
144	IF (HFacW(I,J,K,bi,bj) .NE. 0. D0 ) THEN
145	recip_HFacW(I,J,K,bi,bj) = 1. D0 / HFacW(I,J,K,bi,bj)
146	maskW(I,J,K,bi,bj) = 1. D0
147	ELSE
148	recip_HFacW(I,J,K,bi,bj) = 0. D0
149	maskW(I,J,K,bi,bj) = 0.0 D0
150	ENDIF
151	IF (HFacS(I,J,K,bi,bj) .NE. 0. D0 ) THEN
152	recip_HFacS(I,J,K,bi,bj) = 1. D0 / HFacS(I,J,K,bi,bj)
153	maskS(I,J,K,bi,bj) = 1. D0
154	ELSE
155	recip_HFacS(I,J,K,bi,bj) = 0. D0
156	maskS(I,J,K,bi,bj) = 0. D0
157	ENDIF
158	ENDDO
159	ENDDO
160	ENDDO
161	ENDDO
162	ENDDO
163	_EXCH_XYZ_R4(recip_HFacC , myThid )
164	_EXCH_XYZ_R4(recip_HFacW , myThid )
165	_EXCH_XYZ_R4(recip_HFacS , myThid )
166	_EXCH_XYZ_R4(maskW , myThid )
167	_EXCH_XYZ_R4(maskS , myThid )
168
169	C Calculate recipricols grid lengths
170	DO bj = myByLo(myThid), myByHi(myThid)
171	DO bi = myBxLo(myThid), myBxHi(myThid)
172	DO J=1,sNy
173	DO I=1,sNx
174	recip_dxG(I,J,bi,bj)=1.d0/dxG(I,J,bi,bj)
175	recip_dyG(I,J,bi,bj)=1.d0/dyG(I,J,bi,bj)
176	recip_dxC(I,J,bi,bj)=1.d0/dxC(I,J,bi,bj)
177	recip_dyC(I,J,bi,bj)=1.d0/dyC(I,J,bi,bj)
178	recip_dxF(I,J,bi,bj)=1.d0/dxF(I,J,bi,bj)
179	recip_dyF(I,J,bi,bj)=1.d0/dyF(I,J,bi,bj)
180	recip_dxV(I,J,bi,bj)=1.d0/dxV(I,J,bi,bj)
181	recip_dyU(I,J,bi,bj)=1.d0/dyU(I,J,bi,bj)
182	ENDDO
183	ENDDO
184	ENDDO
185	ENDDO
186	_EXCH_XY_R4(recip_dxG, myThid )
187	_EXCH_XY_R4(recip_dyG, myThid )
188	_EXCH_XY_R4(recip_dxC, myThid )
189	_EXCH_XY_R4(recip_dyC, myThid )
190	_EXCH_XY_R4(recip_dxF, myThid )
191	_EXCH_XY_R4(recip_dyF, myThid )
192	_EXCH_XY_R4(recip_dxV, myThid )
193	_EXCH_XY_R4(recip_dyU, myThid )
194
195	C
196	RETURN
197	END