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

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

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

Revision 1.26 - (show annotations) (download)
Wed Sep 18 16:38:02 2002 UTC (21 years, 7 months ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint46n_post, checkpoint47e_post, checkpoint46l_post, checkpoint47c_post, checkpoint50c_post, checkpoint48e_post, checkpoint50c_pre, checkpoint48i_post, checkpoint46l_pre, checkpoint51, checkpoint50, checkpoint50d_post, checkpoint50b_pre, checkpoint51f_post, checkpoint48b_post, checkpoint51d_post, checkpoint48c_pre, checkpoint47d_pre, checkpoint47a_post, checkpoint48d_pre, checkpoint51j_post, checkpoint47i_post, checkpoint47d_post, checkpoint48d_post, checkpoint48f_post, checkpoint46j_pre, checkpoint48h_post, checkpoint51b_pre, checkpoint47g_post, checkpoint46j_post, checkpoint51h_pre, checkpoint46k_post, checkpoint48a_post, checkpoint50f_post, checkpoint50a_post, checkpoint50f_pre, checkpoint47j_post, branch-exfmods-tag, branchpoint-genmake2, checkpoint48c_post, checkpoint51b_post, checkpoint51c_post, checkpoint47b_post, checkpoint46h_pre, checkpoint46m_post, checkpoint50g_post, checkpoint46g_post, checkpoint50h_post, checkpoint50e_pre, checkpoint50i_post, checkpoint51i_pre, checkpoint47f_post, checkpoint50e_post, checkpoint46i_post, checkpoint50d_pre, checkpoint51e_post, checkpoint47, checkpoint48, checkpoint49, checkpoint46h_post, checkpoint51f_pre, checkpoint48g_post, checkpoint47h_post, checkpoint51g_post, checkpoint50b_post, checkpoint51a_post
Branch point for: branch-exfmods-curt, branch-genmake2
Changes since 1.25: +7 -1 lines 
o Include a new diagnostic variable phiHydLow for the ocean model
  - in z-coordinates, it is the bottom pressure anomaly
  - in p-coordinates, it is the sea surface elevation
  - in both cases, these variable have global drift, reflecting the mass
    drift in z-coordinates and the volume drift in p-coordinates
  - included time averaging for phiHydLow, be aware of the drift!
o depth-dependent computation of Bo_surf for pressure coordinates
  in the ocean (buoyancyRelation='OCEANICP')
  - requires a new routine (FIND_RHO_SCALAR) to compute density with only
    Theta, Salinity, and Pressure in the parameter list. This routine is
    presently contained in find_rho.F. This routine does not give the
    correct density for 'POLY3', which would be a z-dependent reference
    density.
o cleaned up find_rho
  - removed obsolete 'eqn' from the parameter list.
o added two new verification experiments: gop and goz
  (4x4 degree global ocean, 15 layers in pressure and height coordinates)
1 C $Header: /u/gcmpack/MITgcm/model/src/ini_masks_etc.F,v 1.25 2001/11/08 16:36:12 jmc Exp $
2 C $Name: $
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 maskH(i,j,bi,bj) = 0.
132 Ro_surf(I,J,bi,bj) = R_low(I,J,bi,bj)
133 DO K=Nr,1,-1
134 Ro_surf(I,J,bi,bj) = Ro_surf(I,J,bi,bj)
135 & + drF(k)*hFacC(I,J,K,bi,bj)
136 IF (hFacC(I,J,K,bi,bj).NE.0.) THEN
137 ksurfC(I,J,bi,bj) = k
138 maskH(i,j,bi,bj) = 1.
139 tmpfld(i,j,bi,bj) = tmpfld(i,j,bi,bj) + 1.
140 ENDIF
141 ENDDO
142 kLowC(I,J,bi,bj) = 0
143 DO K= 1, Nr
144 IF (hFacC(I,J,K,bi,bj).NE.0) THEN
145 kLowC(I,J,bi,bj) = K
146 ENDIF
147 ENDDO
148 ENDDO
149 ENDDO
150 C - end bi,bj loops.
151 ENDDO
152 ENDDO
153
154 C CALL PLOT_FIELD_XYRS( tmpfld,
155 C & 'Model Depths K Index' , 1, myThid )
156 CALL PLOT_FIELD_XYRS(R_low,
157 & 'Model R_low (ini_masks_etc)', 1, myThid)
158 CALL PLOT_FIELD_XYRS(Ro_surf,
159 & 'Model Ro_surf (ini_masks_etc)', 1, myThid)
160
161 C Calculate quantities derived from XY depth map
162 DO bj = myByLo(myThid), myByHi(myThid)
163 DO bi = myBxLo(myThid), myBxHi(myThid)
164 DO j=1-Oly,sNy+Oly
165 DO i=1-Olx,sNx+Olx
166 C Total fluid column thickness (r_unit) :
167 c Rcolumn(i,j,bi,bj)= Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
168 tmpfld(i,j,bi,bj) = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
169 C Inverse of fluid column thickness (1/r_unit)
170 IF ( tmpfld(i,j,bi,bj) .LE. 0. ) THEN
171 recip_Rcol(i,j,bi,bj) = 0.
172 ELSE
173 recip_Rcol(i,j,bi,bj) = 1. / tmpfld(i,j,bi,bj)
174 ENDIF
175 ENDDO
176 ENDDO
177 ENDDO
178 ENDDO
179 C _EXCH_XY_R4( recip_Rcol, myThid )
180
181 C hFacW and hFacS (at U and V points)
182 DO bj=myByLo(myThid), myByHi(myThid)
183 DO bi=myBxLo(myThid), myBxHi(myThid)
184 DO K=1, Nr
185 DO J=1,sNy
186 DO I=1,sNx
187 hFacW(I,J,K,bi,bj)=
188 & MIN(hFacC(I,J,K,bi,bj),hFacC(I-1,J,K,bi,bj))
189 hFacS(I,J,K,bi,bj)=
190 & MIN(hFacC(I,J,K,bi,bj),hFacC(I,J-1,K,bi,bj))
191 ENDDO
192 ENDDO
193 ENDDO
194 ENDDO
195 ENDDO
196 CALL EXCH_UV_XYZ_RS(hFacW,hFacS,.FALSE.,myThid)
197 C The following block allows thin walls representation of non-periodic
198 C boundaries such as happen on the lat-lon grid at the N/S poles.
199 C We should really supply a flag for doing this.
200 DO bj=myByLo(myThid), myByHi(myThid)
201 DO bi=myBxLo(myThid), myBxHi(myThid)
202 DO K=1, Nr
203 DO J=1-Oly,sNy+Oly
204 DO I=1-Olx,sNx+Olx
205 IF (DYG(I,J,bi,bj).EQ.0.) hFacW(I,J,K,bi,bj)=0.
206 IF (DXG(I,J,bi,bj).EQ.0.) hFacS(I,J,K,bi,bj)=0.
207 ENDDO
208 ENDDO
209 ENDDO
210 ENDDO
211 ENDDO
212
213 C- Write to disk: Total Column Thickness & hFac(C,W,S):
214 _BARRIER
215 _BEGIN_MASTER( myThid )
216 C CALL MDSWRITEFIELD( 'Depth', writeBinaryPrec, .TRUE.,
217 C & 'RS', 1, tmpfld, 1, -1, myThid )
218 CALL WRITE_FLD_XY_RS( 'Depth',' ',tmpfld,0,myThid)
219 CALL WRITE_FLD_XYZ_RS( 'hFacC',' ',hFacC,0,myThid)
220 CALL WRITE_FLD_XYZ_RS( 'hFacW',' ',hFacW,0,myThid)
221 CALL WRITE_FLD_XYZ_RS( 'hFacS',' ',hFacS,0,myThid)
222 _END_MASTER(myThid)
223
224 CALL PLOT_FIELD_XYZRS( hFacC, 'hFacC' , Nr, 1, myThid )
225 CALL PLOT_FIELD_XYZRS( hFacW, 'hFacW' , Nr, 1, myThid )
226 CALL PLOT_FIELD_XYZRS( hFacS, 'hFacS' , Nr, 1, myThid )
227
228 C Masks and reciprocals of hFac[CWS]
229 DO bj = myByLo(myThid), myByHi(myThid)
230 DO bi = myBxLo(myThid), myBxHi(myThid)
231 DO K=1,Nr
232 DO J=1-Oly,sNy+Oly
233 DO I=1-Olx,sNx+Olx
234 IF (HFacC(I,J,K,bi,bj) .NE. 0. ) THEN
235 recip_HFacC(I,J,K,bi,bj) = 1. / HFacC(I,J,K,bi,bj)
236 maskC(I,J,K,bi,bj) = 1.
237 ELSE
238 recip_HFacC(I,J,K,bi,bj) = 0.
239 maskC(I,J,K,bi,bj) = 0.
240 ENDIF
241 IF (HFacW(I,J,K,bi,bj) .NE. 0. ) THEN
242 recip_HFacW(I,J,K,bi,bj) = 1. / HFacW(I,J,K,bi,bj)
243 maskW(I,J,K,bi,bj) = 1.
244 ELSE
245 recip_HFacW(I,J,K,bi,bj) = 0.
246 maskW(I,J,K,bi,bj) = 0.
247 ENDIF
248 IF (HFacS(I,J,K,bi,bj) .NE. 0. ) THEN
249 recip_HFacS(I,J,K,bi,bj) = 1. / HFacS(I,J,K,bi,bj)
250 maskS(I,J,K,bi,bj) = 1.
251 ELSE
252 recip_HFacS(I,J,K,bi,bj) = 0.
253 maskS(I,J,K,bi,bj) = 0.
254 ENDIF
255 ENDDO
256 ENDDO
257 ENDDO
258 C- Calculate surface k index for interface W & S (U & V points)
259 DO J=1-Oly,sNy+Oly
260 DO I=1-Olx,sNx+Olx
261 ksurfW(I,J,bi,bj) = Nr+1
262 ksurfS(I,J,bi,bj) = Nr+1
263 DO k=Nr,1,-1
264 IF (hFacW(I,J,K,bi,bj).NE.0.) ksurfW(I,J,bi,bj) = k
265 IF (hFacS(I,J,K,bi,bj).NE.0.) ksurfS(I,J,bi,bj) = k
266 ENDDO
267 ENDDO
268 ENDDO
269 C - end bi,bj loops.
270 ENDDO
271 ENDDO
272 C _EXCH_XYZ_R4(recip_HFacC , myThid )
273 C _EXCH_XYZ_R4(recip_HFacW , myThid )
274 C _EXCH_XYZ_R4(recip_HFacS , myThid )
275 C _EXCH_XYZ_R4(maskW , myThid )
276 C _EXCH_XYZ_R4(maskS , myThid )
277
278 C Calculate recipricols grid lengths
279 DO bj = myByLo(myThid), myByHi(myThid)
280 DO bi = myBxLo(myThid), myBxHi(myThid)
281 DO J=1-Oly,sNy+Oly
282 DO I=1-Olx,sNx+Olx
283 IF ( dxG(I,J,bi,bj) .NE. 0. )
284 & recip_dxG(I,J,bi,bj)=1.d0/dxG(I,J,bi,bj)
285 IF ( dyG(I,J,bi,bj) .NE. 0. )
286 & recip_dyG(I,J,bi,bj)=1.d0/dyG(I,J,bi,bj)
287 IF ( dxC(I,J,bi,bj) .NE. 0. )
288 & recip_dxC(I,J,bi,bj)=1.d0/dxC(I,J,bi,bj)
289 IF ( dyC(I,J,bi,bj) .NE. 0. )
290 & recip_dyC(I,J,bi,bj)=1.d0/dyC(I,J,bi,bj)
291 IF ( dxF(I,J,bi,bj) .NE. 0. )
292 & recip_dxF(I,J,bi,bj)=1.d0/dxF(I,J,bi,bj)
293 IF ( dyF(I,J,bi,bj) .NE. 0. )
294 & recip_dyF(I,J,bi,bj)=1.d0/dyF(I,J,bi,bj)
295 IF ( dxV(I,J,bi,bj) .NE. 0. )
296 & recip_dxV(I,J,bi,bj)=1.d0/dxV(I,J,bi,bj)
297 IF ( dyU(I,J,bi,bj) .NE. 0. )
298 & recip_dyU(I,J,bi,bj)=1.d0/dyU(I,J,bi,bj)
299 IF ( rA(I,J,bi,bj) .NE. 0. )
300 & recip_rA(I,J,bi,bj)=1.d0/rA(I,J,bi,bj)
301 IF ( rAs(I,J,bi,bj) .NE. 0. )
302 & recip_rAs(I,J,bi,bj)=1.d0/rAs(I,J,bi,bj)
303 IF ( rAw(I,J,bi,bj) .NE. 0. )
304 & recip_rAw(I,J,bi,bj)=1.d0/rAw(I,J,bi,bj)
305 IF ( rAz(I,J,bi,bj) .NE. 0. )
306 & recip_rAz(I,J,bi,bj)=1.d0/rAz(I,J,bi,bj)
307 ENDDO
308 ENDDO
309 ENDDO
310 ENDDO
311 C Do not need these since above denominators are valid over full range
312 C _EXCH_XY_R4(recip_dxG, myThid )
313 C _EXCH_XY_R4(recip_dyG, myThid )
314 C _EXCH_XY_R4(recip_dxC, myThid )
315 C _EXCH_XY_R4(recip_dyC, myThid )
316 C _EXCH_XY_R4(recip_dxF, myThid )
317 C _EXCH_XY_R4(recip_dyF, myThid )
318 C _EXCH_XY_R4(recip_dxV, myThid )
319 C _EXCH_XY_R4(recip_dyU, myThid )
320 C _EXCH_XY_R4(recip_rAw, myThid )
321 C _EXCH_XY_R4(recip_rAs, myThid )
322
323 #ifdef ALLOW_NONHYDROSTATIC
324 C-- Calculate the reciprocal hfac distance/volume for W cells
325 DO bj = myByLo(myThid), myByHi(myThid)
326 DO bi = myBxLo(myThid), myBxHi(myThid)
327 DO K=1,Nr
328 Km1=max(K-1,1)
329 hFacUpper=drF(Km1)/(drF(Km1)+drF(K))
330 IF (Km1.EQ.K) hFacUpper=0.
331 hFacLower=drF(K)/(drF(Km1)+drF(K))
332 DO J=1-Oly,sNy+Oly
333 DO I=1-Olx,sNx+Olx
334 IF (hFacC(I,J,K,bi,bj).NE.0.) THEN
335 IF (hFacC(I,J,K,bi,bj).LE.0.5) THEN
336 recip_hFacU(I,J,K,bi,bj)=
337 & hFacUpper+hFacLower*hFacC(I,J,K,bi,bj)
338 ELSE
339 recip_hFacU(I,J,K,bi,bj)=1.
340 ENDIF
341 ELSE
342 recip_hFacU(I,J,K,bi,bj)=0.
343 ENDIF
344 IF (recip_hFacU(I,J,K,bi,bj).NE.0.)
345 & recip_hFacU(I,J,K,bi,bj)=1./recip_hFacU(I,J,K,bi,bj)
346 ENDDO
347 ENDDO
348 ENDDO
349 ENDDO
350 ENDDO
351 C _EXCH_XY_R4(recip_hFacU, myThid )
352 #endif
353 C
354 RETURN
355 END
  ViewVC Help
Powered by ViewVC 1.1.22