/[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.28 - (hide annotations) (download)
Thu Jan 8 15:28:45 2004 UTC (20 years, 4 months ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint52l_pre, hrcube4, hrcube5, checkpoint52j_pre, checkpoint52l_post, checkpoint52k_post, checkpoint53, checkpoint52f_post, checkpoint52i_pre, hrcube_1, hrcube_2, hrcube_3, checkpoint52m_post, checkpoint52f_pre, checkpoint53a_post, checkpoint52i_post, checkpoint52h_pre, checkpoint52j_post, checkpoint52n_post
Changes since 1.27: +7 -4 lines
Extending loop range for hFacW/S for more generality.
 - This changes nothing and as JMC is keen to point will it not save
   us any problems for a long time but nevertheless it is worth keeping
   improvements you've already made so you don't have to do them again
   in several years time.

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

  ViewVC Help
Powered by ViewVC 1.1.22