/[MITgcm]/MITgcm/model/src/ini_masks_etc.F
ViewVC logotype

Diff of /MITgcm/model/src/ini_masks_etc.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.18 by adcroft, Tue Apr 11 13:39:09 2000 UTC revision 1.19 by adcroft, Fri Feb 2 21:04:48 2001 UTC
# Line 32  C     bi,bj  - Loop counters Line 32  C     bi,bj  - Loop counters
32  C     I,J,K  C     I,J,K
33        INTEGER bi, bj        INTEGER bi, bj
34        INTEGER  I, J, K        INTEGER  I, J, K
35        _RL hFacTmp        INTEGER kadj_rf, klev_noH
36  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
37        INTEGER Km1        INTEGER Km1
38        _RL hFacUpper,hFacLower        _RL hFacUpper,hFacLower
39  #endif  #endif
40          _RL hFacCtmp,topo_rkfac
41          _RL hFacMnSz
42    
43          IF (groundAtK1) THEN
44           topo_rkfac = -rkFac
45           kadj_rf = 1
46           klev_noH = Nr+1
47          ELSE
48           topo_rkfac = rkFac
49           kadj_rf = 0
50           klev_noH = 1
51          ENDIF
52    
53  C     Calculate lopping factor hFacC  C     Calculate lopping factor hFacC
54        DO bj=myByLo(myThid), myByHi(myThid)        DO bj=myByLo(myThid), myByHi(myThid)
55         DO bi=myBxLo(myThid), myBxHi(myThid)         DO bi=myBxLo(myThid), myBxHi(myThid)
56          DO K=1, Nr          DO K=1, Nr
57           DO J=1,sNy           DO J=1-Oly,sNy+Oly
58            DO I=1,sNx            DO I=1-Olx,sNx+Olx
59  C          Round depths within a small fraction of layer depth to that  c          IF (groundAtK1) THEN
60  C          layer depth.  C          o Non-dimensional distance between grid boundary and model depth
61             IF ( ABS(H(I,J,bi,bj)-rF(K)) .LT.  C            for case with "ground" at K=1 (i.e. like a good atmos model)
62       &          1. _d -6*ABS(rF(K)) .AND.  C          e.g. rkfac=+1 => dR/dk<0 (eg. R=p): hFacCtmp=(H(x,y)-rF(k))/drF(K)
63       &          ABS(H(I,J,bi,bj)-rF(K)) .LT.  C          e.g. rkfac=-1 => dR/dk>0 (eg. R=z): hFacCtmp=(rF(K)-H(x,y))/drF(K)
64       &          1. _d -6*ABS(H(I,J,bi,bj)) )THEN  c          hFacCtmp=rkFac*(H(I,J,bi,bj)-rF(K+1))*recip_drF(K)
65              H(I,J,bi,bj) = rF(K)  c          ELSE
66             ENDIF  C          o Non-dimensional distance between grid boundary and model depth
67             IF     ( H(I,J,bi,bj)*rkFac .GE. rF(K)*rkFac ) THEN  C            for case with "ground" at K=Nr (i.e. like original ocean model)
68  C           Top of cell is below base of domain  C          e.g. rkfac=+1 => dR/dk<0 (eg. R=z): hFacCtmp=(rF(K)-H(x,y))/drF(K)
69              hFacC(I,J,K,bi,bj) = 0.  C          e.g. rkfac=-1 => dR/dk>0 (eg. R=p): hFacCtmp=(H(x,y)-rF(k))/drF(K)
70             ELSEIF ( H(I,J,bi,bj)*rkFac .LE. rF(K+1)*rkFac ) THEN  c          hFacCtmp=rkFac*(rF(K)-H(I,J,bi,bj))*recip_drF(K)
71  C           Base of domain is below bottom of this cell  c          ENDIF
72              hFacC(I,J,K,bi,bj) = 1.             hFacCtmp=topo_rkfac*(rF(K+kadj_rf)-H(I,J,bi,bj))*recip_drF(K)
73             ELSE  C          o Select between, closed, open or partial (0,1,0-1)
74  C           Base of domain is in this cell             hFacCtmp=min( max( hFacCtmp, 0. _d 0) , 1. _d 0)
75  C           Set hFac to the fraction of the cell that is open.  C          o And there we have it, the fractional open cell volume
76              hFacC(I,J,K,bi,bj) =             hFacC(I,J,K,bi,bj)=hFacCtmp
77       &       (rF(K)*rkFac-H(I,J,bi,bj)*rkFac)*recip_drF(K)  C          o Impose minimum fraction and/or size (dimensional)
78             ENDIF             hFacMnSz=max( hFacMin , min(hFacMinDr*recip_drF(k),1. _d 0) )
79  C Impose minimum fraction and/or size             IF (hFacC(I,J,K,bi,bj).LT.hFacMnSz) THEN
80             hFacTmp=max( hFacMin , min(hFacMinDr*recip_drF(k),1.) )              IF (hFacC(I,J,K,bi,bj).LT.hFacMnSz*0.5) THEN
            IF (hFacC(I,J,K,bi,bj).LT.hFacTmp) THEN  
             IF (hFacC(I,J,K,bi,bj).LT.hFacTmp*0.5) THEN  
81               hFacC(I,J,K,bi,bj)=0.               hFacC(I,J,K,bi,bj)=0.
82              ELSE              ELSE
83               hFacC(I,J,K,bi,bj)=hFacTmp               hFacC(I,J,K,bi,bj)=hFacMnSz
84              ENDIF              ENDIF
85             ENDIF             ENDIF
86             depthInK(i,j,bi,bj) = depthInK(i,j,bi,bj) + 1.             IF (hFacC(I,J,K,bi,bj).NE.0.)
87  Crg  &                          +hFacC(i,j,k,bi,bj)       &         depthInK(i,j,bi,bj) = depthInK(i,j,bi,bj) + 1.
88            ENDDO            ENDDO
89           ENDDO           ENDDO
90          ENDDO          ENDDO
91         ENDDO         ENDDO
92        ENDDO        ENDDO
93        _EXCH_XYZ_R4(hFacC , myThid )  C     _EXCH_XYZ_R4(hFacC , myThid )
94        _EXCH_XY_R4( depthInK, myThid )  C     _EXCH_XY_R4( depthInK, myThid )
95    
96        CALL PLOT_FIELD_XYRS( depthInK,        CALL PLOT_FIELD_XYRS( depthInK,
97       & 'Model Depths K Index' , 1, myThid )       & 'Model Depths K Index' , 1, myThid )
# Line 89  Crg  &                          +hFacC(i Line 99  Crg  &                          +hFacC(i
99  C Re-calculate depth of ocean, taking into account hFacC  C Re-calculate depth of ocean, taking into account hFacC
100        DO bj=myByLo(myThid), myByHi(myThid)        DO bj=myByLo(myThid), myByHi(myThid)
101         DO bi=myBxLo(myThid), myBxHi(myThid)         DO bi=myBxLo(myThid), myBxHi(myThid)
102          DO J=1,sNy          DO J=1-Oly,sNy+Oly
103           DO I=1,sNx           DO I=1-Olx,sNx+Olx
104            H(I,J,bi,bj)=0.            H(I,J,bi,bj)=rF(klev_noH)
105            DO K=1,Nr            DO K=1,Nr
106             H(I,J,bi,bj)=H(I,J,bi,bj)-             H(I,J,bi,bj)=H(I,J,bi,bj)-
107       &       rkFac*drF(k)*hFacC(I,J,K,bi,bj)       &       topo_rkFac*drF(k)*hFacC(I,J,K,bi,bj)
108            ENDDO            ENDDO
109           ENDDO           ENDDO
110          ENDDO          ENDDO
111         ENDDO         ENDDO
112        ENDDO        ENDDO
113        _EXCH_XY_R4(H , myThid )  C     _EXCH_XY_R4(H , myThid )
114          CALL PLOT_FIELD_XYRS(H,'Model depths (ini_masks_etc)',1,myThid)
115        CALL WRITE_FLD_XY_RS( 'Depth',' ',H,0,myThid)        CALL WRITE_FLD_XY_RS( 'Depth',' ',H,0,myThid)
116          CALL WRITE_FLD_XYZ_RS( 'hFacC',' ',hFacC,0,myThid)
117  C     CALL MDSWRITEFIELD( 'Depth', writeBinaryPrec, .TRUE.,  C     CALL MDSWRITEFIELD( 'Depth', writeBinaryPrec, .TRUE.,
118  C    &                    'RS', 1, H, 1, -1, myThid )  C    &                    'RS', 1, H, 1, -1, myThid )
119    
120  C     Calculate quantities derived from XY depth map  C     Calculate quantities derived from XY depth map
121        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
122         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
123          DO J=1,sNy          DO J=1-Oly,sNy+Oly
124           DO I=1,sNx           DO I=1-Olx,sNx+Olx
125  C         Inverse of depth  C         Inverse of depth
126            IF ( h(i,j,bi,bj) .EQ. 0. _d 0 ) THEN            IF ( h(i,j,bi,bj) .EQ. 0. ) THEN
127             recip_H(i,j,bi,bj) = 0. _d 0             recip_H(i,j,bi,bj) = 0.
128            ELSE            ELSE
129             recip_H(i,j,bi,bj) = 1. _d 0 /  abs( H(i,j,bi,bj) )             recip_H(i,j,bi,bj) = 1. / abs( H(i,j,bi,bj) )
130            ENDIF            ENDIF
131            depthInK(i,j,bi,bj) = 0.            depthInK(i,j,bi,bj) = 0.
132           ENDDO           ENDDO
133          ENDDO          ENDDO
134         ENDDO         ENDDO
135        ENDDO        ENDDO
136        _EXCH_XY_R4(   recip_H, myThid )  C     _EXCH_XY_R4(   recip_H, myThid )
137    
138  C     hFacW and hFacS (at U and V points)  C     hFacW and hFacS (at U and V points)
139        DO bj=myByLo(myThid), myByHi(myThid)        DO bj=myByLo(myThid), myByHi(myThid)
# Line 140  C     hFacW and hFacS (at U and V points Line 152  C     hFacW and hFacS (at U and V points
152        ENDDO        ENDDO
153        _EXCH_XYZ_R4(hFacW , myThid )        _EXCH_XYZ_R4(hFacW , myThid )
154        _EXCH_XYZ_R4(hFacS , myThid )        _EXCH_XYZ_R4(hFacS , myThid )
155    C     Re-do hFacW and hFacS (at U and V points)
156          DO bj=myByLo(myThid), myByHi(myThid)
157           DO bi=myBxLo(myThid), myBxHi(myThid)
158            DO K=1, Nr
159             DO J=1-Oly,sNy+Oly
160              DO I=1-Olx+1,sNx+Olx ! Note range
161               hFacW(I,J,K,bi,bj)=
162         &       MIN(hFacC(I,J,K,bi,bj),hFacC(I-1,J,K,bi,bj))
163              ENDDO
164             ENDDO
165             DO I=1-Olx,sNx+Olx
166              DO J=1-Oly+1,sNy+Oly ! Note range
167               hFacS(I,J,K,bi,bj)=
168         &       MIN(hFacC(I,J,K,bi,bj),hFacC(I,J-1,K,bi,bj))
169              ENDDO
170             ENDDO
171            ENDDO
172           ENDDO
173          ENDDO
174    
175          CALL PLOT_FIELD_XYZRS( hFacC, 'hFacC' , Nr, 1, myThid )
176          CALL PLOT_FIELD_XYZRS( hFacW, 'hFacW' , Nr, 1, myThid )
177          CALL PLOT_FIELD_XYZRS( hFacS, 'hFacS' , Nr, 1, myThid )
178    
179  C     Masks and reciprocals of hFac[CWS]  C     Masks and reciprocals of hFac[CWS]
180        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
181         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
182          DO K=1,Nr          DO K=1,Nr
183           DO J=1,sNy           DO J=1-Oly,sNy+Oly
184            DO I=1,sNx            DO I=1-Olx,sNx+Olx
185             IF (HFacC(I,J,K,bi,bj) .NE. 0. _d 0 ) THEN             IF (HFacC(I,J,K,bi,bj) .NE. 0. ) THEN
186              recip_HFacC(I,J,K,bi,bj) = 1. _d 0 / HFacC(I,J,K,bi,bj)              recip_HFacC(I,J,K,bi,bj) = 1. / HFacC(I,J,K,bi,bj)
187             ELSE             ELSE
188              recip_HFacC(I,J,K,bi,bj) = 0. _d 0              recip_HFacC(I,J,K,bi,bj) = 0.
189             ENDIF             ENDIF
190             IF (HFacW(I,J,K,bi,bj) .NE. 0. _d 0 ) THEN             IF (HFacW(I,J,K,bi,bj) .NE. 0. ) THEN
191              recip_HFacW(I,J,K,bi,bj) = 1. _d 0 / HFacW(I,J,K,bi,bj)              recip_HFacW(I,J,K,bi,bj) = 1. / HFacW(I,J,K,bi,bj)
192              maskW(I,J,K,bi,bj) = 1. _d 0              maskW(I,J,K,bi,bj) = 1.
193             ELSE             ELSE
194              recip_HFacW(I,J,K,bi,bj) = 0. _d 0              recip_HFacW(I,J,K,bi,bj) = 0.
195              maskW(I,J,K,bi,bj) = 0.0 _d 0              maskW(I,J,K,bi,bj) = 0.
196             ENDIF             ENDIF
197             IF (HFacS(I,J,K,bi,bj) .NE. 0. _d 0 ) THEN             IF (HFacS(I,J,K,bi,bj) .NE. 0. ) THEN
198              recip_HFacS(I,J,K,bi,bj) = 1. _d 0 / HFacS(I,J,K,bi,bj)              recip_HFacS(I,J,K,bi,bj) = 1. / HFacS(I,J,K,bi,bj)
199              maskS(I,J,K,bi,bj) = 1. _d 0              maskS(I,J,K,bi,bj) = 1.
200             ELSE             ELSE
201              recip_HFacS(I,J,K,bi,bj) = 0. _d 0              recip_HFacS(I,J,K,bi,bj) = 0.
202              maskS(I,J,K,bi,bj) = 0. _d 0              maskS(I,J,K,bi,bj) = 0.
203             ENDIF             ENDIF
204            ENDDO            ENDDO
205           ENDDO           ENDDO
206          ENDDO          ENDDO
207         ENDDO         ENDDO
208        ENDDO        ENDDO
209        _EXCH_XYZ_R4(recip_HFacC    , myThid )  C     _EXCH_XYZ_R4(recip_HFacC    , myThid )
210        _EXCH_XYZ_R4(recip_HFacW    , myThid )  C     _EXCH_XYZ_R4(recip_HFacW    , myThid )
211        _EXCH_XYZ_R4(recip_HFacS    , myThid )  C     _EXCH_XYZ_R4(recip_HFacS    , myThid )
212        _EXCH_XYZ_R4(maskW    , myThid )  C     _EXCH_XYZ_R4(maskW    , myThid )
213        _EXCH_XYZ_R4(maskS    , myThid )  C     _EXCH_XYZ_R4(maskS    , myThid )
214    
215  C     Calculate recipricols grid lengths  C     Calculate recipricols grid lengths
216        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
217         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
218          DO J=1,sNy          DO J=1-Oly,sNy+Oly
219           DO I=1,sNx           DO I=1-Olx,sNx+Olx
220            recip_dxG(I,J,bi,bj)=1.d0/dxG(I,J,bi,bj)            IF ( dxG(I,J,bi,bj) .NE. 0. )
221            recip_dyG(I,J,bi,bj)=1.d0/dyG(I,J,bi,bj)       &    recip_dxG(I,J,bi,bj)=1.d0/dxG(I,J,bi,bj)
222            recip_dxC(I,J,bi,bj)=1.d0/dxC(I,J,bi,bj)            IF ( dyG(I,J,bi,bj) .NE. 0. )
223            recip_dyC(I,J,bi,bj)=1.d0/dyC(I,J,bi,bj)       &    recip_dyG(I,J,bi,bj)=1.d0/dyG(I,J,bi,bj)
224            recip_dxF(I,J,bi,bj)=1.d0/dxF(I,J,bi,bj)            IF ( dxC(I,J,bi,bj) .NE. 0. )
225            recip_dyF(I,J,bi,bj)=1.d0/dyF(I,J,bi,bj)       &    recip_dxC(I,J,bi,bj)=1.d0/dxC(I,J,bi,bj)
226            recip_dxV(I,J,bi,bj)=1.d0/dxV(I,J,bi,bj)            IF ( dyC(I,J,bi,bj) .NE. 0. )
227            recip_dyU(I,J,bi,bj)=1.d0/dyU(I,J,bi,bj)       &    recip_dyC(I,J,bi,bj)=1.d0/dyC(I,J,bi,bj)
228              IF ( dxF(I,J,bi,bj) .NE. 0. )
229         &    recip_dxF(I,J,bi,bj)=1.d0/dxF(I,J,bi,bj)
230              IF ( dyF(I,J,bi,bj) .NE. 0. )
231         &    recip_dyF(I,J,bi,bj)=1.d0/dyF(I,J,bi,bj)
232              IF ( dxV(I,J,bi,bj) .NE. 0. )
233         &    recip_dxV(I,J,bi,bj)=1.d0/dxV(I,J,bi,bj)
234              IF ( dyU(I,J,bi,bj) .NE. 0. )
235         &    recip_dyU(I,J,bi,bj)=1.d0/dyU(I,J,bi,bj)
236              IF ( rA(I,J,bi,bj) .NE. 0. )
237         &    recip_rA(I,J,bi,bj)=1.d0/rA(I,J,bi,bj)
238              IF ( rAs(I,J,bi,bj) .NE. 0. )
239         &    recip_rAs(I,J,bi,bj)=1.d0/rAs(I,J,bi,bj)
240              IF ( rAw(I,J,bi,bj) .NE. 0. )
241         &    recip_rAw(I,J,bi,bj)=1.d0/rAw(I,J,bi,bj)
242              IF ( rAz(I,J,bi,bj) .NE. 0. )
243         &    recip_rAz(I,J,bi,bj)=1.d0/rAz(I,J,bi,bj)
244           ENDDO           ENDDO
245          ENDDO          ENDDO
246         ENDDO         ENDDO
247        ENDDO        ENDDO
248        _EXCH_XY_R4(recip_dxG, myThid )  C     Do not need these since above denominators are valid over full range
249        _EXCH_XY_R4(recip_dyG, myThid )  C     _EXCH_XY_R4(recip_dxG, myThid )
250        _EXCH_XY_R4(recip_dxC, myThid )  C     _EXCH_XY_R4(recip_dyG, myThid )
251        _EXCH_XY_R4(recip_dyC, myThid )  C     _EXCH_XY_R4(recip_dxC, myThid )
252        _EXCH_XY_R4(recip_dxF, myThid )  C     _EXCH_XY_R4(recip_dyC, myThid )
253        _EXCH_XY_R4(recip_dyF, myThid )  C     _EXCH_XY_R4(recip_dxF, myThid )
254        _EXCH_XY_R4(recip_dxV, myThid )  C     _EXCH_XY_R4(recip_dyF, myThid )
255        _EXCH_XY_R4(recip_dyU, myThid )  C     _EXCH_XY_R4(recip_dxV, myThid )
256    C     _EXCH_XY_R4(recip_dyU, myThid )
257    C     _EXCH_XY_R4(recip_rAw, myThid )
258    C     _EXCH_XY_R4(recip_rAs, myThid )
259    
260  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
261  C--   Calculate the reciprocal hfac distance/volume for W cells  C--   Calculate the reciprocal hfac distance/volume for W cells
# Line 212  C--   Calculate the reciprocal hfac dist Line 266  C--   Calculate the reciprocal hfac dist
266           hFacUpper=drF(Km1)/(drF(Km1)+drF(K))           hFacUpper=drF(Km1)/(drF(Km1)+drF(K))
267           IF (Km1.EQ.K) hFacUpper=0.           IF (Km1.EQ.K) hFacUpper=0.
268           hFacLower=drF(K)/(drF(Km1)+drF(K))           hFacLower=drF(K)/(drF(Km1)+drF(K))
269           DO J=1,sNy           DO J=1-Oly,sNy+Oly
270            DO I=1,sNx            DO I=1-Olx,sNx+Olx
271             IF (hFacC(I,J,K,bi,bj).NE.0.) THEN             IF (hFacC(I,J,K,bi,bj).NE.0.) THEN
272              IF (hFacC(I,J,K,bi,bj).LE.0.5) THEN              IF (hFacC(I,J,K,bi,bj).LE.0.5) THEN
273              recip_hFacU(I,J,K,bi,bj)=              recip_hFacU(I,J,K,bi,bj)=
# Line 231  C--   Calculate the reciprocal hfac dist Line 285  C--   Calculate the reciprocal hfac dist
285          ENDDO          ENDDO
286         ENDDO         ENDDO
287        ENDDO        ENDDO
288        _EXCH_XY_R4(recip_hFacU, myThid )  C     _EXCH_XY_R4(recip_hFacU, myThid )
289  #endif  #endif
290  C  C
291        RETURN        RETURN
Legend:
Removed from v.1.18  
changed lines
  Added in v.1.19
  ViewVC Help
Powered by ViewVC 1.1.22