/[MITgcm]/MITgcm_contrib/AITCZ/code/ini_masks_etc.F
ViewVC logotype

Annotation of /MITgcm_contrib/AITCZ/code/ini_masks_etc.F

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


Revision 1.1 - (hide annotations) (download)
Wed Aug 20 15:24:59 2003 UTC (21 years, 11 months ago) by czaja
Branch: MAIN
CVS Tags: HEAD
Initial creation of Arnaud's simple coupled simulation.

1 czaja 1.1 C $Header: /u/u0/gcmpack/models/MITgcmUV/model/src/ini_masks_etc.F,v 1.24 2001/09/26 18:09:15 cnh Exp $
2     C $Name: release1_beta1 $
3    
4     #include "CPP_OPTIONS.h"
5    
6     CBOP
7     C !ROUTINE: INI_MASKS_ETC
8     C !INTERFACE:
9     SUBROUTINE INI_MASKS_ETC( myThid )
10     C !DESCRIPTION: \bv
11     C *==========================================================*
12     C | SUBROUTINE INI_MASKS_ETC
13     C | o Initialise masks and topography factors
14     C *==========================================================*
15     C | These arrays are used throughout the code and describe
16     C | the topography of the domain through masks (0s and 1s)
17     C | and fractional height factors (0<hFac<1). The latter
18     C | distinguish between the lopped-cell and full-step
19     C | topographic representations.
20     C *==========================================================*
21     C \ev
22    
23     C !USES:
24     IMPLICIT NONE
25     C === Global variables ===
26     #include "SIZE.h"
27     #include "EEPARAMS.h"
28     #include "PARAMS.h"
29     #include "GRID.h"
30     #include "SURFACE.h"
31    
32     C !INPUT/OUTPUT PARAMETERS:
33     C == Routine arguments ==
34     C myThid - Number of this instance of INI_MASKS_ETC
35     INTEGER myThid
36    
37     C !LOCAL VARIABLES:
38     C == Local variables in common ==
39     C tmpfld - Temporary array used to compute & write Total Depth
40     C has to be in common for multi threading
41     COMMON / LOCAL_INI_MASKS_ETC / tmpfld
42     _RS tmpfld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
43     C == Local variables ==
44     C bi,bj - Loop counters
45     C I,J,K
46     INTEGER bi, bj
47     INTEGER I, J, K
48     #ifdef ALLOW_NONHYDROSTATIC
49     INTEGER Km1
50     _RL hFacUpper,hFacLower
51     #endif
52     _RL hFacCtmp
53     _RL hFacMnSz
54     CEOP
55    
56     C- Calculate lopping factor hFacC : over-estimate the part inside of the domain
57     C taking into account the lower_R Boundary (Bathymetrie / Top of Atmos)
58     DO bj=myByLo(myThid), myByHi(myThid)
59     DO bi=myBxLo(myThid), myBxHi(myThid)
60     DO K=1, Nr
61     hFacMnSz=max( hFacMin, min(hFacMinDr*recip_drF(k),1. _d 0) )
62     DO J=1-Oly,sNy+Oly
63     DO I=1-Olx,sNx+Olx
64     C o Non-dimensional distance between grid bound. and domain lower_R bound.
65     hFacCtmp = (rF(K)-R_low(I,J,bi,bj))*recip_drF(K)
66     C o Select between, closed, open or partial (0,1,0-1)
67     hFacCtmp=min( max( hFacCtmp, 0. _d 0) , 1. _d 0)
68     C o Impose minimum fraction and/or size (dimensional)
69     IF (hFacCtmp.LT.hFacMnSz) THEN
70     IF (hFacCtmp.LT.hFacMnSz*0.5) THEN
71     hFacC(I,J,K,bi,bj)=0.
72     ELSE
73     hFacC(I,J,K,bi,bj)=hFacMnSz
74     ENDIF
75     ELSE
76     hFacC(I,J,K,bi,bj)=hFacCtmp
77     ENDIF
78     ENDDO
79     ENDDO
80     ENDDO
81    
82     C- Re-calculate lower-R Boundary position, taking into account hFacC
83     DO J=1-Oly,sNy+Oly
84     DO I=1-Olx,sNx+Olx
85     R_low(I,J,bi,bj) = rF(1)
86     DO K=Nr,1,-1
87     R_low(I,J,bi,bj) = R_low(I,J,bi,bj)
88     & - drF(k)*hFacC(I,J,K,bi,bj)
89     ENDDO
90     ENDDO
91     ENDDO
92     C - end bi,bj loops.
93     ENDDO
94     ENDDO
95    
96     C- Calculate lopping factor hFacC : Remove part outside of the domain
97     C taking into account the Reference (=at rest) Surface Position Ro_surf
98     DO bj=myByLo(myThid), myByHi(myThid)
99     DO bi=myBxLo(myThid), myBxHi(myThid)
100     DO K=1, Nr
101     hFacMnSz=max( hFacMin, min(hFacMinDr*recip_drF(k),1. _d 0) )
102     DO J=1-Oly,sNy+Oly
103     DO I=1-Olx,sNx+Olx
104     C o Non-dimensional distance between grid boundary and model surface
105     hFacCtmp = (rF(k)-Ro_surf(I,J,bi,bj))*recip_drF(K)
106     C o Reduce the previous fraction : substract the outside part.
107     hFacCtmp = hFacC(I,J,K,bi,bj) - max( hFacCtmp, 0. _d 0)
108     C o set to zero if empty Column :
109     hFacCtmp = max( hFacCtmp, 0. _d 0)
110     C o Impose minimum fraction and/or size (dimensional)
111     IF (hFacCtmp.LT.hFacMnSz) THEN
112     IF (hFacCtmp.LT.hFacMnSz*0.5) THEN
113     hFacC(I,J,K,bi,bj)=0.
114     ELSE
115     hFacC(I,J,K,bi,bj)=hFacMnSz
116     ENDIF
117     ELSE
118     hFacC(I,J,K,bi,bj)=hFacCtmp
119     ENDIF
120     ENDDO
121     ENDDO
122     ENDDO
123    
124     C- Re-calculate Reference surface position, taking into account hFacC
125     C initialize Total column fluid thickness and surface k index
126     C Note: if no fluid (continent) ==> ksurf = Nr+1
127     DO J=1-Oly,sNy+Oly
128     DO I=1-Olx,sNx+Olx
129     tmpfld(I,J,bi,bj) = 0.
130     ksurfC(I,J,bi,bj) = Nr+1
131     Ro_surf(I,J,bi,bj) = R_low(I,J,bi,bj)
132     DO K=Nr,1,-1
133     Ro_surf(I,J,bi,bj) = Ro_surf(I,J,bi,bj)
134     & + drF(k)*hFacC(I,J,K,bi,bj)
135     IF (hFacC(I,J,K,bi,bj).NE.0.) THEN
136     ksurfC(I,J,bi,bj) = k
137     tmpfld(i,j,bi,bj) = tmpfld(i,j,bi,bj) + 1.
138     ENDIF
139     ENDDO
140     ENDDO
141     ENDDO
142     C - end bi,bj loops.
143     ENDDO
144     ENDDO
145    
146     C CALL PLOT_FIELD_XYRS( tmpfld,
147     C & 'Model Depths K Index' , 1, myThid )
148     CALL PLOT_FIELD_XYRS(R_low,
149     & 'Model R_low (ini_masks_etc)', 1, myThid)
150     CALL PLOT_FIELD_XYRS(Ro_surf,
151     & 'Model Ro_surf (ini_masks_etc)', 1, myThid)
152    
153     C Calculate quantities derived from XY depth map
154     DO bj = myByLo(myThid), myByHi(myThid)
155     DO bi = myBxLo(myThid), myBxHi(myThid)
156     DO j=1-Oly,sNy+Oly
157     DO i=1-Olx,sNx+Olx
158     C Total fluid column thickness (r_unit) :
159     c Rcolumn(i,j,bi,bj)= Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
160     tmpfld(i,j,bi,bj) = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
161     C Inverse of fluid column thickness (1/r_unit)
162     IF ( tmpfld(i,j,bi,bj) .LE. 0. ) THEN
163     recip_Rcol(i,j,bi,bj) = 0.
164     ELSE
165     recip_Rcol(i,j,bi,bj) = 1. / tmpfld(i,j,bi,bj)
166     ENDIF
167     ENDDO
168     ENDDO
169     ENDDO
170     ENDDO
171     C _EXCH_XY_R4( recip_Rcol, myThid )
172    
173     C hFacW and hFacS (at U and V points)
174     DO bj=myByLo(myThid), myByHi(myThid)
175     DO bi=myBxLo(myThid), myBxHi(myThid)
176     DO K=1, Nr
177     DO J=1,sNy
178     DO I=1,sNx
179     hFacW(I,J,K,bi,bj)=
180     & MIN(hFacC(I,J,K,bi,bj),hFacC(I-1,J,K,bi,bj))
181     hFacS(I,J,K,bi,bj)=
182     & MIN(hFacC(I,J,K,bi,bj),hFacC(I,J-1,K,bi,bj))
183     ENDDO
184     ENDDO
185     ENDDO
186     ENDDO
187     ENDDO
188     CALL EXCH_UV_XYZ_RS(hFacW,hFacS,.FALSE.,myThid)
189     C The following block allows thin walls representation of non-periodic
190     C boundaries such as happen on the lat-lon grid at the N/S poles.
191     C We should really supply a flag for doing this.
192     DO bj=myByLo(myThid), myByHi(myThid)
193     DO bi=myBxLo(myThid), myBxHi(myThid)
194     DO K=1, Nr
195     DO J=1-Oly,sNy+Oly
196     DO I=1-Olx,sNx+Olx
197     IF (DYG(I,J,bi,bj).EQ.0.) hFacW(I,J,K,bi,bj)=0.
198     IF (DXG(I,J,bi,bj).EQ.0.) hFacS(I,J,K,bi,bj)=0.
199     ENDDO
200     ENDDO
201     ENDDO
202     ENDDO
203     ENDDO
204    
205     C- Write to disk: Total Column Thickness & hFac(C,W,S):
206     _BARRIER
207     _BEGIN_MASTER( myThid )
208     C CALL MDSWRITEFIELD( 'Depth', writeBinaryPrec, .TRUE.,
209     C & 'RS', 1, tmpfld, 1, -1, myThid )
210     CC(acz) CALL WRITE_FLD_XY_RS( 'Depth',' ',tmpfld,0,myThid)
211     CC(acz) CALL WRITE_FLD_XYZ_RS( 'hFacC',' ',hFacC,0,myThid)
212     CC(acz) CALL WRITE_FLD_XYZ_RS( 'hFacW',' ',hFacW,0,myThid)
213     CC(acz) CALL WRITE_FLD_XYZ_RS( 'hFacS',' ',hFacS,0,myThid)
214     _END_MASTER(myThid)
215    
216     CALL PLOT_FIELD_XYZRS( hFacC, 'hFacC' , Nr, 1, myThid )
217     CALL PLOT_FIELD_XYZRS( hFacW, 'hFacW' , Nr, 1, myThid )
218     CALL PLOT_FIELD_XYZRS( hFacS, 'hFacS' , Nr, 1, myThid )
219    
220     C Masks and reciprocals of hFac[CWS]
221     DO bj = myByLo(myThid), myByHi(myThid)
222     DO bi = myBxLo(myThid), myBxHi(myThid)
223     DO K=1,Nr
224     DO J=1-Oly,sNy+Oly
225     DO I=1-Olx,sNx+Olx
226     IF (HFacC(I,J,K,bi,bj) .NE. 0. ) THEN
227     recip_HFacC(I,J,K,bi,bj) = 1. / HFacC(I,J,K,bi,bj)
228     maskC(I,J,K,bi,bj) = 1.
229     ELSE
230     recip_HFacC(I,J,K,bi,bj) = 0.
231     maskC(I,J,K,bi,bj) = 0.
232     ENDIF
233     IF (HFacW(I,J,K,bi,bj) .NE. 0. ) THEN
234     recip_HFacW(I,J,K,bi,bj) = 1. / HFacW(I,J,K,bi,bj)
235     maskW(I,J,K,bi,bj) = 1.
236     ELSE
237     recip_HFacW(I,J,K,bi,bj) = 0.
238     maskW(I,J,K,bi,bj) = 0.
239     ENDIF
240     IF (HFacS(I,J,K,bi,bj) .NE. 0. ) THEN
241     recip_HFacS(I,J,K,bi,bj) = 1. / HFacS(I,J,K,bi,bj)
242     maskS(I,J,K,bi,bj) = 1.
243     ELSE
244     recip_HFacS(I,J,K,bi,bj) = 0.
245     maskS(I,J,K,bi,bj) = 0.
246     ENDIF
247     ENDDO
248     ENDDO
249     ENDDO
250     C- Calculate surface k index for interface W & S (U & V points)
251     DO J=1-Oly,sNy+Oly
252     DO I=1-Olx,sNx+Olx
253     ksurfW(I,J,bi,bj) = Nr+1
254     ksurfS(I,J,bi,bj) = Nr+1
255     DO k=Nr,1,-1
256     IF (hFacW(I,J,K,bi,bj).NE.0.) ksurfW(I,J,bi,bj) = k
257     IF (hFacS(I,J,K,bi,bj).NE.0.) ksurfS(I,J,bi,bj) = k
258     ENDDO
259     ENDDO
260     ENDDO
261     C - end bi,bj loops.
262     ENDDO
263     ENDDO
264     C _EXCH_XYZ_R4(recip_HFacC , myThid )
265     C _EXCH_XYZ_R4(recip_HFacW , myThid )
266     C _EXCH_XYZ_R4(recip_HFacS , myThid )
267     C _EXCH_XYZ_R4(maskW , myThid )
268     C _EXCH_XYZ_R4(maskS , myThid )
269    
270     C Calculate recipricols grid lengths
271     DO bj = myByLo(myThid), myByHi(myThid)
272     DO bi = myBxLo(myThid), myBxHi(myThid)
273     DO J=1-Oly,sNy+Oly
274     DO I=1-Olx,sNx+Olx
275     IF ( dxG(I,J,bi,bj) .NE. 0. )
276     & recip_dxG(I,J,bi,bj)=1.d0/dxG(I,J,bi,bj)
277     IF ( dyG(I,J,bi,bj) .NE. 0. )
278     & recip_dyG(I,J,bi,bj)=1.d0/dyG(I,J,bi,bj)
279     IF ( dxC(I,J,bi,bj) .NE. 0. )
280     & recip_dxC(I,J,bi,bj)=1.d0/dxC(I,J,bi,bj)
281     IF ( dyC(I,J,bi,bj) .NE. 0. )
282     & recip_dyC(I,J,bi,bj)=1.d0/dyC(I,J,bi,bj)
283     IF ( dxF(I,J,bi,bj) .NE. 0. )
284     & recip_dxF(I,J,bi,bj)=1.d0/dxF(I,J,bi,bj)
285     IF ( dyF(I,J,bi,bj) .NE. 0. )
286     & recip_dyF(I,J,bi,bj)=1.d0/dyF(I,J,bi,bj)
287     IF ( dxV(I,J,bi,bj) .NE. 0. )
288     & recip_dxV(I,J,bi,bj)=1.d0/dxV(I,J,bi,bj)
289     IF ( dyU(I,J,bi,bj) .NE. 0. )
290     & recip_dyU(I,J,bi,bj)=1.d0/dyU(I,J,bi,bj)
291     IF ( rA(I,J,bi,bj) .NE. 0. )
292     & recip_rA(I,J,bi,bj)=1.d0/rA(I,J,bi,bj)
293     IF ( rAs(I,J,bi,bj) .NE. 0. )
294     & recip_rAs(I,J,bi,bj)=1.d0/rAs(I,J,bi,bj)
295     IF ( rAw(I,J,bi,bj) .NE. 0. )
296     & recip_rAw(I,J,bi,bj)=1.d0/rAw(I,J,bi,bj)
297     IF ( rAz(I,J,bi,bj) .NE. 0. )
298     & recip_rAz(I,J,bi,bj)=1.d0/rAz(I,J,bi,bj)
299     ENDDO
300     ENDDO
301     ENDDO
302     ENDDO
303     C Do not need these since above denominators are valid over full range
304     C _EXCH_XY_R4(recip_dxG, myThid )
305     C _EXCH_XY_R4(recip_dyG, myThid )
306     C _EXCH_XY_R4(recip_dxC, myThid )
307     C _EXCH_XY_R4(recip_dyC, myThid )
308     C _EXCH_XY_R4(recip_dxF, myThid )
309     C _EXCH_XY_R4(recip_dyF, myThid )
310     C _EXCH_XY_R4(recip_dxV, myThid )
311     C _EXCH_XY_R4(recip_dyU, myThid )
312     C _EXCH_XY_R4(recip_rAw, myThid )
313     C _EXCH_XY_R4(recip_rAs, myThid )
314    
315     #ifdef ALLOW_NONHYDROSTATIC
316     C-- Calculate the reciprocal hfac distance/volume for W cells
317     DO bj = myByLo(myThid), myByHi(myThid)
318     DO bi = myBxLo(myThid), myBxHi(myThid)
319     DO K=1,Nr
320     Km1=max(K-1,1)
321     hFacUpper=drF(Km1)/(drF(Km1)+drF(K))
322     IF (Km1.EQ.K) hFacUpper=0.
323     hFacLower=drF(K)/(drF(Km1)+drF(K))
324     DO J=1-Oly,sNy+Oly
325     DO I=1-Olx,sNx+Olx
326     IF (hFacC(I,J,K,bi,bj).NE.0.) THEN
327     IF (hFacC(I,J,K,bi,bj).LE.0.5) THEN
328     recip_hFacU(I,J,K,bi,bj)=
329     & hFacUpper+hFacLower*hFacC(I,J,K,bi,bj)
330     ELSE
331     recip_hFacU(I,J,K,bi,bj)=1.
332     ENDIF
333     ELSE
334     recip_hFacU(I,J,K,bi,bj)=0.
335     ENDIF
336     IF (recip_hFacU(I,J,K,bi,bj).NE.0.)
337     & recip_hFacU(I,J,K,bi,bj)=1./recip_hFacU(I,J,K,bi,bj)
338     ENDDO
339     ENDDO
340     ENDDO
341     ENDDO
342     ENDDO
343     C _EXCH_XY_R4(recip_hFacU, myThid )
344     #endif
345     C
346     RETURN
347     END

  ViewVC Help
Powered by ViewVC 1.1.22