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

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

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


Revision 1.31 - (hide annotations) (download)
Tue Feb 7 11:47:48 2006 UTC (18 years, 4 months ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint58b_post, checkpoint58f_post, checkpoint58d_post, checkpoint58a_post, checkpoint58m_post, checkpoint58e_post, checkpoint58k_post, checkpoint58l_post, checkpoint58g_post, checkpoint58h_post, checkpoint58j_post, checkpoint58i_post, checkpoint58c_post
Changes since 1.30: +33 -1 lines
o add hooks for new package shelfice, painless

1 mlosch 1.31 C $Header: /u/gcmpack/MITgcm/model/src/ini_masks_etc.F,v 1.30 2005/07/20 22:31:47 jmc 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     C | SUBROUTINE INI_MASKS_ETC
14     C | o Initialise masks and topography factors
15     C *==========================================================*
16     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     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 mlosch 1.31 #ifdef ALLOW_SHELFICE
33     # include "SHELFICE.h"
34     #endif /* ALLOW_SHELFICE */
35 adcroft 1.1
36 cnh 1.24 C !INPUT/OUTPUT PARAMETERS:
37 adcroft 1.1 C == Routine arguments ==
38 cnh 1.6 C myThid - Number of this instance of INI_MASKS_ETC
39 adcroft 1.1 INTEGER myThid
40    
41 cnh 1.24 C !LOCAL VARIABLES:
42 jmc 1.22 C == Local variables in common ==
43     C tmpfld - Temporary array used to compute & write Total Depth
44     C has to be in common for multi threading
45     COMMON / LOCAL_INI_MASKS_ETC / tmpfld
46     _RS tmpfld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
47 adcroft 1.1 C == Local variables ==
48     C bi,bj - Loop counters
49     C I,J,K
50     INTEGER bi, bj
51     INTEGER I, J, K
52 adcroft 1.15 #ifdef ALLOW_NONHYDROSTATIC
53     INTEGER Km1
54     _RL hFacUpper,hFacLower
55     #endif
56 adcroft 1.21 _RL hFacCtmp
57 adcroft 1.19 _RL hFacMnSz
58 jmc 1.30 _RL tileArea
59 cnh 1.24 CEOP
60 adcroft 1.19
61 adcroft 1.21 C- Calculate lopping factor hFacC : over-estimate the part inside of the domain
62     C taking into account the lower_R Boundary (Bathymetrie / Top of Atmos)
63 adcroft 1.2 DO bj=myByLo(myThid), myByHi(myThid)
64     DO bi=myBxLo(myThid), myBxHi(myThid)
65 cnh 1.4 DO K=1, Nr
66 adcroft 1.21 hFacMnSz=max( hFacMin, min(hFacMinDr*recip_drF(k),1. _d 0) )
67 adcroft 1.19 DO J=1-Oly,sNy+Oly
68     DO I=1-Olx,sNx+Olx
69 adcroft 1.21 C o Non-dimensional distance between grid bound. and domain lower_R bound.
70     hFacCtmp = (rF(K)-R_low(I,J,bi,bj))*recip_drF(K)
71     C o Select between, closed, open or partial (0,1,0-1)
72 adcroft 1.19 hFacCtmp=min( max( hFacCtmp, 0. _d 0) , 1. _d 0)
73 adcroft 1.21 C o Impose minimum fraction and/or size (dimensional)
74     IF (hFacCtmp.LT.hFacMnSz) THEN
75     IF (hFacCtmp.LT.hFacMnSz*0.5) THEN
76 adcroft 1.3 hFacC(I,J,K,bi,bj)=0.
77     ELSE
78 adcroft 1.19 hFacC(I,J,K,bi,bj)=hFacMnSz
79 adcroft 1.3 ENDIF
80 adcroft 1.21 ELSE
81     hFacC(I,J,K,bi,bj)=hFacCtmp
82 adcroft 1.2 ENDIF
83     ENDDO
84     ENDDO
85     ENDDO
86 adcroft 1.21
87     C- Re-calculate lower-R Boundary position, taking into account hFacC
88     DO J=1-Oly,sNy+Oly
89     DO I=1-Olx,sNx+Olx
90     R_low(I,J,bi,bj) = rF(1)
91     DO K=Nr,1,-1
92     R_low(I,J,bi,bj) = R_low(I,J,bi,bj)
93     & - drF(k)*hFacC(I,J,K,bi,bj)
94     ENDDO
95     ENDDO
96     ENDDO
97     C - end bi,bj loops.
98 adcroft 1.2 ENDDO
99     ENDDO
100 cnh 1.7
101 adcroft 1.21 C- Calculate lopping factor hFacC : Remove part outside of the domain
102     C taking into account the Reference (=at rest) Surface Position Ro_surf
103 adcroft 1.16 DO bj=myByLo(myThid), myByHi(myThid)
104     DO bi=myBxLo(myThid), myBxHi(myThid)
105 adcroft 1.21 DO K=1, Nr
106     hFacMnSz=max( hFacMin, min(hFacMinDr*recip_drF(k),1. _d 0) )
107     DO J=1-Oly,sNy+Oly
108     DO I=1-Olx,sNx+Olx
109     C o Non-dimensional distance between grid boundary and model surface
110     hFacCtmp = (rF(k)-Ro_surf(I,J,bi,bj))*recip_drF(K)
111     C o Reduce the previous fraction : substract the outside part.
112     hFacCtmp = hFacC(I,J,K,bi,bj) - max( hFacCtmp, 0. _d 0)
113     C o set to zero if empty Column :
114     hFacCtmp = max( hFacCtmp, 0. _d 0)
115     C o Impose minimum fraction and/or size (dimensional)
116     IF (hFacCtmp.LT.hFacMnSz) THEN
117     IF (hFacCtmp.LT.hFacMnSz*0.5) THEN
118     hFacC(I,J,K,bi,bj)=0.
119     ELSE
120     hFacC(I,J,K,bi,bj)=hFacMnSz
121     ENDIF
122     ELSE
123     hFacC(I,J,K,bi,bj)=hFacCtmp
124     ENDIF
125     ENDDO
126     ENDDO
127     ENDDO
128    
129 mlosch 1.31 #ifdef ALLOW_SHELFICE
130     C-- compute contributions of shelf ice to looping factors
131     IF ( useShelfIce ) THEN
132     DO K=1, Nr
133     hFacMnSz=max( hFacMin, min(hFacMinDr*recip_drF(k),1. _d 0) )
134     DO J=1-Oly,sNy+Oly
135     DO I=1-Olx,sNx+Olx
136     C o Non-dimensional distance between grid boundary and model surface
137     hFacCtmp = (rF(k)-R_shelfIce(I,J,bi,bj))*recip_drF(K)
138     C o Reduce the previous fraction : substract the outside part.
139     hFacCtmp = hFacC(I,J,K,bi,bj) - max( hFacCtmp, 0. _d 0)
140     C o set to zero if empty Column :
141     hFacCtmp = max( hFacCtmp, 0. _d 0)
142     C o Impose minimum fraction and/or size (dimensional)
143     IF (hFacCtmp.LT.hFacMnSz) THEN
144     IF (hFacCtmp.LT.hFacMnSz*0.5) THEN
145     hFacC(I,J,K,bi,bj)=0.
146     ELSE
147     hFacC(I,J,K,bi,bj)=hFacMnSz
148     ENDIF
149     ELSE
150     hFacC(I,J,K,bi,bj)=hFacCtmp
151     ENDIF
152     ENDDO
153     ENDDO
154     ENDDO
155     ENDIF
156     #endif /* ALLOW_SHELFICE */
157    
158 adcroft 1.21 C- Re-calculate Reference surface position, taking into account hFacC
159     C initialize Total column fluid thickness and surface k index
160 jmc 1.23 C Note: if no fluid (continent) ==> ksurf = Nr+1
161 adcroft 1.19 DO J=1-Oly,sNy+Oly
162     DO I=1-Olx,sNx+Olx
163 adcroft 1.21 tmpfld(I,J,bi,bj) = 0.
164 jmc 1.23 ksurfC(I,J,bi,bj) = Nr+1
165 jmc 1.25 maskH(i,j,bi,bj) = 0.
166 adcroft 1.21 Ro_surf(I,J,bi,bj) = R_low(I,J,bi,bj)
167     DO K=Nr,1,-1
168     Ro_surf(I,J,bi,bj) = Ro_surf(I,J,bi,bj)
169     & + drF(k)*hFacC(I,J,K,bi,bj)
170     IF (hFacC(I,J,K,bi,bj).NE.0.) THEN
171 jmc 1.23 ksurfC(I,J,bi,bj) = k
172 jmc 1.25 maskH(i,j,bi,bj) = 1.
173 adcroft 1.21 tmpfld(i,j,bi,bj) = tmpfld(i,j,bi,bj) + 1.
174 mlosch 1.26 ENDIF
175     ENDDO
176     kLowC(I,J,bi,bj) = 0
177     DO K= 1, Nr
178     IF (hFacC(I,J,K,bi,bj).NE.0) THEN
179     kLowC(I,J,bi,bj) = K
180 adcroft 1.21 ENDIF
181 adcroft 1.16 ENDDO
182     ENDDO
183     ENDDO
184 adcroft 1.21 C - end bi,bj loops.
185 adcroft 1.16 ENDDO
186     ENDDO
187 adcroft 1.21
188     C CALL PLOT_FIELD_XYRS( tmpfld,
189     C & 'Model Depths K Index' , 1, myThid )
190     CALL PLOT_FIELD_XYRS(R_low,
191     & 'Model R_low (ini_masks_etc)', 1, myThid)
192     CALL PLOT_FIELD_XYRS(Ro_surf,
193     & 'Model Ro_surf (ini_masks_etc)', 1, myThid)
194 adcroft 1.16
195     C Calculate quantities derived from XY depth map
196 jmc 1.30 globalArea = 0. _d 0
197 adcroft 1.16 DO bj = myByLo(myThid), myByHi(myThid)
198     DO bi = myBxLo(myThid), myBxHi(myThid)
199 adcroft 1.21 DO j=1-Oly,sNy+Oly
200     DO i=1-Olx,sNx+Olx
201     C Total fluid column thickness (r_unit) :
202     c Rcolumn(i,j,bi,bj)= Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
203     tmpfld(i,j,bi,bj) = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
204     C Inverse of fluid column thickness (1/r_unit)
205     IF ( tmpfld(i,j,bi,bj) .LE. 0. ) THEN
206     recip_Rcol(i,j,bi,bj) = 0.
207 adcroft 1.16 ELSE
208 adcroft 1.21 recip_Rcol(i,j,bi,bj) = 1. / tmpfld(i,j,bi,bj)
209 adcroft 1.16 ENDIF
210     ENDDO
211     ENDDO
212 jmc 1.30 C- Compute the domain Area:
213     tileArea = 0. _d 0
214     DO j=1,sNy
215     DO i=1,sNx
216     tileArea = tileArea + rA(i,j,bi,bj)*maskH(i,j,bi,bj)
217     ENDDO
218     ENDDO
219     globalArea = globalArea + tileArea
220 adcroft 1.16 ENDDO
221     ENDDO
222 adcroft 1.21 C _EXCH_XY_R4( recip_Rcol, myThid )
223 jmc 1.30 _GLOBAL_SUM_R8( globalArea, myThid )
224 adcroft 1.1
225     C hFacW and hFacS (at U and V points)
226     DO bj=myByLo(myThid), myByHi(myThid)
227     DO bi=myBxLo(myThid), myBxHi(myThid)
228 cnh 1.4 DO K=1, Nr
229 adcroft 1.28 DO J=1-Oly,sNy+Oly
230     DO I=2-Olx,sNx+Olx
231 adcroft 1.1 hFacW(I,J,K,bi,bj)=
232     & MIN(hFacC(I,J,K,bi,bj),hFacC(I-1,J,K,bi,bj))
233 adcroft 1.28 ENDDO
234     ENDDO
235     DO J=2-Oly,sNy+oly
236     DO I=1-Olx,sNx+Olx
237 adcroft 1.1 hFacS(I,J,K,bi,bj)=
238     & MIN(hFacC(I,J,K,bi,bj),hFacC(I,J-1,K,bi,bj))
239     ENDDO
240     ENDDO
241     ENDDO
242     ENDDO
243     ENDDO
244 adcroft 1.21 CALL EXCH_UV_XYZ_RS(hFacW,hFacS,.FALSE.,myThid)
245     C The following block allows thin walls representation of non-periodic
246     C boundaries such as happen on the lat-lon grid at the N/S poles.
247     C We should really supply a flag for doing this.
248 adcroft 1.19 DO bj=myByLo(myThid), myByHi(myThid)
249     DO bi=myBxLo(myThid), myBxHi(myThid)
250     DO K=1, Nr
251     DO J=1-Oly,sNy+Oly
252 adcroft 1.21 DO I=1-Olx,sNx+Olx
253     IF (DYG(I,J,bi,bj).EQ.0.) hFacW(I,J,K,bi,bj)=0.
254     IF (DXG(I,J,bi,bj).EQ.0.) hFacS(I,J,K,bi,bj)=0.
255 adcroft 1.19 ENDDO
256     ENDDO
257     ENDDO
258     ENDDO
259     ENDDO
260 jmc 1.22
261     C- Write to disk: Total Column Thickness & hFac(C,W,S):
262     _BARRIER
263 adcroft 1.29 c _BEGIN_MASTER( myThid )
264     C This I/O is now done in write_grid.F
265 jmc 1.22 C CALL MDSWRITEFIELD( 'Depth', writeBinaryPrec, .TRUE.,
266     C & 'RS', 1, tmpfld, 1, -1, myThid )
267 adcroft 1.29 c CALL WRITE_FLD_XY_RS( 'Depth',' ',tmpfld,0,myThid)
268     c CALL WRITE_FLD_XYZ_RS( 'hFacC',' ',hFacC,0,myThid)
269     c CALL WRITE_FLD_XYZ_RS( 'hFacW',' ',hFacW,0,myThid)
270     c CALL WRITE_FLD_XYZ_RS( 'hFacS',' ',hFacS,0,myThid)
271     c _END_MASTER(myThid)
272 adcroft 1.19
273     CALL PLOT_FIELD_XYZRS( hFacC, 'hFacC' , Nr, 1, myThid )
274     CALL PLOT_FIELD_XYZRS( hFacW, 'hFacW' , Nr, 1, myThid )
275     CALL PLOT_FIELD_XYZRS( hFacS, 'hFacS' , Nr, 1, myThid )
276 adcroft 1.1
277     C Masks and reciprocals of hFac[CWS]
278     DO bj = myByLo(myThid), myByHi(myThid)
279     DO bi = myBxLo(myThid), myBxHi(myThid)
280 cnh 1.4 DO K=1,Nr
281 adcroft 1.19 DO J=1-Oly,sNy+Oly
282     DO I=1-Olx,sNx+Olx
283     IF (HFacC(I,J,K,bi,bj) .NE. 0. ) THEN
284     recip_HFacC(I,J,K,bi,bj) = 1. / HFacC(I,J,K,bi,bj)
285 adcroft 1.21 maskC(I,J,K,bi,bj) = 1.
286 adcroft 1.1 ELSE
287 adcroft 1.19 recip_HFacC(I,J,K,bi,bj) = 0.
288 adcroft 1.21 maskC(I,J,K,bi,bj) = 0.
289 adcroft 1.1 ENDIF
290 adcroft 1.19 IF (HFacW(I,J,K,bi,bj) .NE. 0. ) THEN
291     recip_HFacW(I,J,K,bi,bj) = 1. / HFacW(I,J,K,bi,bj)
292     maskW(I,J,K,bi,bj) = 1.
293 adcroft 1.1 ELSE
294 adcroft 1.19 recip_HFacW(I,J,K,bi,bj) = 0.
295     maskW(I,J,K,bi,bj) = 0.
296 adcroft 1.1 ENDIF
297 adcroft 1.19 IF (HFacS(I,J,K,bi,bj) .NE. 0. ) THEN
298     recip_HFacS(I,J,K,bi,bj) = 1. / HFacS(I,J,K,bi,bj)
299     maskS(I,J,K,bi,bj) = 1.
300 adcroft 1.1 ELSE
301 adcroft 1.19 recip_HFacS(I,J,K,bi,bj) = 0.
302     maskS(I,J,K,bi,bj) = 0.
303 adcroft 1.1 ENDIF
304     ENDDO
305     ENDDO
306     ENDDO
307 jmc 1.23 C- Calculate surface k index for interface W & S (U & V points)
308     DO J=1-Oly,sNy+Oly
309     DO I=1-Olx,sNx+Olx
310     ksurfW(I,J,bi,bj) = Nr+1
311     ksurfS(I,J,bi,bj) = Nr+1
312     DO k=Nr,1,-1
313     IF (hFacW(I,J,K,bi,bj).NE.0.) ksurfW(I,J,bi,bj) = k
314     IF (hFacS(I,J,K,bi,bj).NE.0.) ksurfS(I,J,bi,bj) = k
315     ENDDO
316     ENDDO
317     ENDDO
318     C - end bi,bj loops.
319 adcroft 1.1 ENDDO
320     ENDDO
321 adcroft 1.19 C _EXCH_XYZ_R4(recip_HFacC , myThid )
322     C _EXCH_XYZ_R4(recip_HFacW , myThid )
323     C _EXCH_XYZ_R4(recip_HFacS , myThid )
324     C _EXCH_XYZ_R4(maskW , myThid )
325     C _EXCH_XYZ_R4(maskS , myThid )
326 adcroft 1.1
327     C Calculate recipricols grid lengths
328     DO bj = myByLo(myThid), myByHi(myThid)
329     DO bi = myBxLo(myThid), myBxHi(myThid)
330 adcroft 1.19 DO J=1-Oly,sNy+Oly
331     DO I=1-Olx,sNx+Olx
332     IF ( dxG(I,J,bi,bj) .NE. 0. )
333     & recip_dxG(I,J,bi,bj)=1.d0/dxG(I,J,bi,bj)
334     IF ( dyG(I,J,bi,bj) .NE. 0. )
335     & recip_dyG(I,J,bi,bj)=1.d0/dyG(I,J,bi,bj)
336     IF ( dxC(I,J,bi,bj) .NE. 0. )
337     & recip_dxC(I,J,bi,bj)=1.d0/dxC(I,J,bi,bj)
338     IF ( dyC(I,J,bi,bj) .NE. 0. )
339     & recip_dyC(I,J,bi,bj)=1.d0/dyC(I,J,bi,bj)
340     IF ( dxF(I,J,bi,bj) .NE. 0. )
341     & recip_dxF(I,J,bi,bj)=1.d0/dxF(I,J,bi,bj)
342     IF ( dyF(I,J,bi,bj) .NE. 0. )
343     & recip_dyF(I,J,bi,bj)=1.d0/dyF(I,J,bi,bj)
344     IF ( dxV(I,J,bi,bj) .NE. 0. )
345     & recip_dxV(I,J,bi,bj)=1.d0/dxV(I,J,bi,bj)
346     IF ( dyU(I,J,bi,bj) .NE. 0. )
347     & recip_dyU(I,J,bi,bj)=1.d0/dyU(I,J,bi,bj)
348     IF ( rA(I,J,bi,bj) .NE. 0. )
349     & recip_rA(I,J,bi,bj)=1.d0/rA(I,J,bi,bj)
350     IF ( rAs(I,J,bi,bj) .NE. 0. )
351     & recip_rAs(I,J,bi,bj)=1.d0/rAs(I,J,bi,bj)
352     IF ( rAw(I,J,bi,bj) .NE. 0. )
353     & recip_rAw(I,J,bi,bj)=1.d0/rAw(I,J,bi,bj)
354     IF ( rAz(I,J,bi,bj) .NE. 0. )
355     & recip_rAz(I,J,bi,bj)=1.d0/rAz(I,J,bi,bj)
356 adcroft 1.1 ENDDO
357     ENDDO
358     ENDDO
359     ENDDO
360 adcroft 1.19 C Do not need these since above denominators are valid over full range
361     C _EXCH_XY_R4(recip_dxG, myThid )
362     C _EXCH_XY_R4(recip_dyG, myThid )
363     C _EXCH_XY_R4(recip_dxC, myThid )
364     C _EXCH_XY_R4(recip_dyC, myThid )
365     C _EXCH_XY_R4(recip_dxF, myThid )
366     C _EXCH_XY_R4(recip_dyF, myThid )
367     C _EXCH_XY_R4(recip_dxV, myThid )
368     C _EXCH_XY_R4(recip_dyU, myThid )
369     C _EXCH_XY_R4(recip_rAw, myThid )
370     C _EXCH_XY_R4(recip_rAs, myThid )
371 adcroft 1.1
372 adcroft 1.15 #ifdef ALLOW_NONHYDROSTATIC
373     C-- Calculate the reciprocal hfac distance/volume for W cells
374     DO bj = myByLo(myThid), myByHi(myThid)
375     DO bi = myBxLo(myThid), myBxHi(myThid)
376     DO K=1,Nr
377     Km1=max(K-1,1)
378     hFacUpper=drF(Km1)/(drF(Km1)+drF(K))
379     IF (Km1.EQ.K) hFacUpper=0.
380     hFacLower=drF(K)/(drF(Km1)+drF(K))
381 adcroft 1.19 DO J=1-Oly,sNy+Oly
382     DO I=1-Olx,sNx+Olx
383 adcroft 1.15 IF (hFacC(I,J,K,bi,bj).NE.0.) THEN
384     IF (hFacC(I,J,K,bi,bj).LE.0.5) THEN
385     recip_hFacU(I,J,K,bi,bj)=
386     & hFacUpper+hFacLower*hFacC(I,J,K,bi,bj)
387     ELSE
388     recip_hFacU(I,J,K,bi,bj)=1.
389     ENDIF
390     ELSE
391     recip_hFacU(I,J,K,bi,bj)=0.
392     ENDIF
393     IF (recip_hFacU(I,J,K,bi,bj).NE.0.)
394     & recip_hFacU(I,J,K,bi,bj)=1./recip_hFacU(I,J,K,bi,bj)
395     ENDDO
396     ENDDO
397     ENDDO
398     ENDDO
399     ENDDO
400 adcroft 1.19 C _EXCH_XY_R4(recip_hFacU, myThid )
401 adcroft 1.15 #endif
402 adcroft 1.1 C
403     RETURN
404     END

  ViewVC Help
Powered by ViewVC 1.1.22