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

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

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


Revision 1.22 - (show annotations) (download)
Fri Jul 26 14:22:22 2013 UTC (10 years, 10 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65t, checkpoint65u, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint65l, checkpoint65m, checkpoint65b, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e, checkpoint65, checkpoint65o, checkpoint64y, checkpoint64x, checkpoint64z, checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64w, checkpoint64v, checkpoint64m, checkpoint64l, checkpoint64o, checkpoint64n
Changes since 1.21: +6 -3 lines
extend length of drC & recip_drC from Nr to Nr+1.

1 C $Header: /u/gcmpack/MITgcm/model/src/ini_vertical_grid.F,v 1.21 2011/08/18 20:29:49 jmc Exp $
2 C $Name: $
3
4 #include "CPP_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: INI_VERTICAL_GRID
8 C !INTERFACE:
9 SUBROUTINE INI_VERTICAL_GRID( myThid )
10
11 C !DESCRIPTION: \bv
12 C *==========================================================*
13 C | SUBROUTINE INI_VERTICAL_GRID
14 C | o Initialise vertical gridding arrays
15 C *==========================================================*
16 C \ev
17
18 C !USES:
19 IMPLICIT NONE
20 C === Global variables ===
21 #include "SIZE.h"
22 #include "EEPARAMS.h"
23 #include "PARAMS.h"
24 #include "GRID.h"
25
26 C !INPUT/OUTPUT PARAMETERS:
27 C == Routine arguments ==
28 C myThid :: my Thread Id number
29 INTEGER myThid
30
31 C !LOCAL VARIABLES:
32 C == Local variables ==
33 C k :: loop index
34 C msgBuf :: Informational/error message buffer
35 INTEGER k
36 _RL tmpRatio, checkRatio1, checkRatio2
37 CHARACTER*(MAX_LEN_MBUF) msgBuf
38 _RL maxErrC, maxErrF, epsil, tmpError
39 _RL rFullDepth, recip_fullDepth
40 _RS rSigBndRS, tmpRS
41 CEOP
42
43 _BEGIN_MASTER(myThid)
44
45 WRITE(msgBuf,'(A,2(A,L5))') 'Enter INI_VERTICAL_GRID:',
46 & ' setInterFDr=', setInterFDr,
47 & ' ; setCenterDr=', setCenterDr
48 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
49 & SQUEEZE_RIGHT, myThid )
50
51 C-- Set factors required for mixing pressure and meters as vertical coordinate.
52 C rkSign is a "sign" parameter which is used where the orientation of the vertical
53 C coordinate (pressure or meters) relative to the vertical index (k) is important.
54 C rkSign = -1 applies when k and the coordinate are in the opposite sense.
55 C rkSign = 1 applies when k and the coordinate are in the same sense.
56 rkSign = -1. _d 0
57 gravitySign = -1. _d 0
58 IF ( usingPCoords ) THEN
59 gravitySign = 1. _d 0
60 ENDIF
61
62 IF ( .NOT.(setInterFDr.OR.setCenterDr) ) THEN
63 WRITE(msgBuf,'(A)')
64 & 'S/R INI_VERTICAL_GRID: neither delR nor delRc are defined'
65 CALL PRINT_ERROR( msgBuf, myThid )
66 WRITE(msgBuf,'(A)')
67 & 'S/R INI_VERTICAL_GRID: Need at least 1 of the 2 (delR,delRc)'
68 CALL PRINT_ERROR( msgBuf, myThid )
69 STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID'
70 ENDIF
71
72 C--- Set Level r-thickness (drF) and Center r-distances (drC)
73
74 IF (setInterFDr) THEN
75 C-- Interface r-distances are defined:
76 DO k=1,Nr
77 drF(k) = delR(k)
78 ENDDO
79 C- Check that all thickness are > 0 :
80 DO k=1,Nr
81 IF (delR(k).LE.0.) THEN
82 WRITE(msgBuf,'(A,I4,A,E16.8)')
83 & 'S/R INI_VERTICAL_GRID: delR(k=',k,' )=',delR(k)
84 CALL PRINT_ERROR( msgBuf, myThid )
85 WRITE(msgBuf,'(A)')
86 & 'S/R INI_VERTICAL_GRID: Vert. grid spacing MUST BE > 0'
87 CALL PRINT_ERROR( msgBuf, myThid )
88 STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID'
89 ENDIF
90 ENDDO
91 ELSE
92 C-- Interface r-distances undefined:
93 C assume Interface at middle between 2 Center
94 drF(1) = delRc(1)
95 DO k=2,Nr
96 c drF(k-1) = 0.5 _d 0 *delRc(k) + drF(k-1)
97 c drF( k ) = 0.5 _d 0 *delRc(k)
98 C- note: change the order to prevent some compilers to produce wrong code
99 C when trying to optimise this loop :
100 drF( k ) = 0.5 _d 0 *delRc(k)
101 drF(k-1) = 0.5 _d 0 *delRc(k) + drF(k-1)
102 ENDDO
103 drF(Nr) = delRc(Nr+1) + drF(Nr)
104 ENDIF
105
106 IF (setCenterDr) THEN
107 C-- Cell Center r-distances are defined:
108 DO k=1,Nr+1
109 drC(k) = delRc(k)
110 ENDDO
111 C- Check that all thickness are > 0 :
112 DO k=1,Nr+1
113 IF (delRc(k).LE.0.) THEN
114 WRITE(msgBuf,'(A,I4,A,E16.8)')
115 & 'S/R INI_VERTICAL_GRID: delRc(k=',k,' )=',delRc(k)
116 CALL PRINT_ERROR( msgBuf, myThid )
117 WRITE(msgBuf,'(A)')
118 & 'S/R INI_VERTICAL_GRID: Vert. grid spacing MUST BE > 0'
119 CALL PRINT_ERROR( msgBuf, myThid )
120 STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID'
121 ENDIF
122 ENDDO
123 ELSE
124 C-- Cell Center r-distances undefined:
125 C assume Center at middle between 2 Interfaces
126 drC(1) = 0.5 _d 0 *delR(1)
127 DO k=2,Nr
128 drC(k) = 0.5 _d 0 *(delR(k-1)+delR(k))
129 ENDDO
130 drC(Nr+1) = 0.5 _d 0 *delR(Nr)
131 ENDIF
132
133 C--- Set r-position of interFace (rF) and cell-Center (rC):
134 rF(1) = Ro_SeaLevel
135 DO k=1,Nr
136 rF(k+1) = rF(k) + rkSign*drF(k)
137 ENDDO
138 rC(1) = rF(1) + rkSign*drC(1)
139 DO k=2,Nr
140 rC(k) = rC(k-1) + rkSign*drC(k)
141 ENDDO
142
143 C--- Check vertical discretization :
144 checkRatio2 = 100.
145 checkRatio1 = 1. _d 0 / checkRatio2
146 DO k=1,Nr
147 tmpRatio = 0.
148 IF ( (rC(k)-rF(k+1)) .NE. 0. )
149 & tmpRatio = (rF(k)-rC(k)) / (rC(k)-rF(k+1))
150 IF ( tmpRatio.LT.checkRatio1 .OR. tmpRatio.GT.checkRatio2 ) THEN
151 c write(0,*) 'drF=',drF
152 c write(0,*) 'drC=',drC
153 c write(0,*) 'rF=',rF
154 c write(0,*) 'rC=',rC
155 WRITE(msgBuf,'(A,I4,A,E16.8)')
156 & 'S/R INI_VERTICAL_GRID: Invalid relative position, level k=',
157 & k, ' :', tmpRatio
158 CALL PRINT_ERROR( msgBuf, myThid )
159 WRITE(msgBuf,'(A,1PE14.6,A,2E14.6)')
160 & 'S/R INI_VERTICAL_GRID: rC=', rC(k),
161 & ' , rF(k,k+1)=',rF(k),rF(k+1)
162 CALL PRINT_ERROR( msgBuf, myThid )
163 STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID'
164 ENDIF
165 ENDDO
166
167 C- Calculate reciprol vertical grid spacing :
168 DO k=1,Nr+1
169 recip_drC(k) = 1. _d 0/drC(k)
170 ENDDO
171 DO k=1,Nr
172 recip_drF(k) = 1. _d 0/drF(k)
173 ENDDO
174
175 C--- Hybrid-Sigma vertical coordinate:
176 IF ( selectSigmaCoord .EQ. 0 ) THEN
177 DO k=1,Nr+1
178 aHybSigmF(k) = 0. _d 0
179 bHybSigmF(k) = 0. _d 0
180 dAHybSigC(k) = 0. _d 0
181 dAHybSigC(k) = 0. _d 0
182 ENDDO
183 DO k=1,Nr
184 aHybSigmC(k) = 0. _d 0
185 bHybSigmC(k) = 0. _d 0
186 dAHybSigF(k) = 0. _d 0
187 dAHybSigF(k) = 0. _d 0
188 ENDDO
189 ELSE
190 rFullDepth = rF(1) - rF(Nr+1)
191 recip_fullDepth = 0. _d 0
192 IF ( rFullDepth.GT.0. ) recip_fullDepth = 1. _d 0 / rFullDepth
193 rSigBndRS = rSigmaBnd
194 IF ( hybSigmFile.EQ.' ' .AND. rSigmaBnd.EQ.UNSET_RL ) THEN
195 C- Default is pure sigma:
196 IF ( usingPCoords ) rSigBndRS = rF(Nr+1)
197 IF ( usingZCoords ) rSigBndRS = rF(1)
198 ENDIF
199 c IF ( hybSigmFile.EQ.' ' .AND. rSigmaBnd.EQ.UNSET_RL ) THEN
200 C- compute coeff as pure sigma coordinate
201 c DO k=1,Nr+1
202 c aHybSigmF(k) = 0.
203 c bHybSigmF(k) = (rF(k)-rF(Nr+1))*recip_fullDepth
204 c ENDDO
205 c DO k=1,Nr
206 c aHybSigmC(k) = 0.
207 c bHybSigmC(k) = (rC(k)-rF(Nr+1))*recip_fullDepth
208 c ENDDO
209 c ELSEIF ( hybSigmFile.EQ.' ' ) THEN
210 IF ( hybSigmFile.EQ.' ' ) THEN
211 C-- compute coeff assuming fixed r-coord above rSigmaBnd and pure sigma below
212 IF ( usingPCoords .AND. setInterFDr ) THEN
213 C- Alternative method : p = pTop + A*DeltaFullP + B*(eta+Pr_surf - rSigmaBnd)
214 c aHybSigmF(k) = ( MIN(rF(k),rSigmaBnd) - rF(Nr+1) )
215 c & *recip_fullDepth
216 c bHybSigmF(k) = ( MAX(rF(k),rSigmaBnd) - rSigmaBnd )
217 c & /( rF(1) - rSigmaBnd )
218 C- Standard method : p = pTop + A*DeltaFullP + B*(eta+Ro_surf - pTop)
219 C sigma part goes from 0 @ rSigmaBnd (and above) to 1 @ surface
220 DO k=1,Nr+1
221 C- separate the 2 cases:
222 c IF ( rF(k).LE.rSigmaBnd ) THEN
223 c bHybSigmF(k) = 0.
224 c aHybSigmF(k) = ( rF(k) - rF(Nr+1) )*recip_fullDepth
225 c ELSE
226 c bHybSigmF(k) = (rF(k)-rSigmaBnd)/(rF(1)-rSigmaBnd)
227 c aHybSigmF(k) = (1. _d 0 - bHybSigmF(k))
228 c & *(rSigmaBnd-rF(Nr+1) )*recip_fullDepth
229 c ENDIF
230 C- unique formula using min fct:
231 tmpRS = MIN( rF(k), rSigBndRS )
232 bHybSigmF(k) = ( rF(k) - tmpRS )/(rF(1)-rSigBndRS)
233 aHybSigmF(k) = (1. _d 0 - bHybSigmF(k))
234 & *( tmpRS -rF(Nr+1) )*recip_fullDepth
235 ENDDO
236 ENDIF
237 IF ( usingPCoords .AND. setCenterDr ) THEN
238 DO k=1,Nr
239 C- separate the 2 cases:
240 c IF ( rC(k).LE.rSigmaBnd ) THEN
241 c bHybSigmC(k) = 0.
242 c aHybSigmC(k) = ( rC(k) - rF(Nr+1) )*recip_fullDepth
243 c ELSE
244 c bHybSigmC(k) = (rC(k)-rSigmaBnd)/(rF(1)-rSigmaBnd)
245 c aHybSigmC(k) = (1. _d 0 - bHybSigmC(k))
246 c & *(rSigmaBnd-rF(Nr+1) )*recip_fullDepth
247 c ENDIF
248 C- unique formula using min fct:
249 tmpRS = MIN( rC(k), rSigBndRS )
250 bHybSigmC(k) = ( rC(k) - tmpRS )/(rF(1)-rSigBndRS)
251 aHybSigmC(k) = (1. _d 0 - bHybSigmC(k))
252 & *( tmpRS -rF(Nr+1) )*recip_fullDepth
253 ENDDO
254 ENDIF
255 IF ( usingZCoords .AND. setInterFDr ) THEN
256 C- Standard method : z = zBot + A*DeltaFullZ + B*(eta+Ro_surf - zBot)
257 C sigma part goes from 1 @ rSigmaBnd (and above) to 0 @ bottom
258 DO k=1,Nr+1
259 C- separate the 2 cases:
260 c IF ( rF(k).GE.rSigmaBnd ) THEN
261 c bHybSigmF(k) = 1.
262 c aHybSigmF(k) = ( rF(k)-rF(1) )*recip_fullDepth
263 c ELSE
264 c bHybSigmF(k) = ( rF(k)-rF(Nr+1) )/( rSigmaBnd-rF(Nr+1) )
265 c aHybSigmF(k) = bHybSigmF(k)*(rSigmaBnd-rF(1))*recip_fullDepth
266 c ENDIF
267 C- unique formula using max fct:
268 tmpRS = MAX( rF(k), rSigBndRS )
269 bHybSigmF(k) = ( rF(k)-rF(Nr+1) )/( tmpRS-rF(Nr+1) )
270 aHybSigmF(k) = bHybSigmF(k)*( tmpRS-rF(1) )*recip_fullDepth
271 ENDDO
272 ENDIF
273 IF ( usingZCoords .AND. setCenterDr ) THEN
274 DO k=1,Nr
275 C- separate the 2 cases:
276 c IF ( rC(k).GE.rSigmaBnd ) THEN
277 c bHybSigmC(k) = 1.
278 c aHybSigmC(k) = ( rC(k)-rF(1) )*recip_fullDepth
279 c ELSE
280 c bHybSigmC(k) = ( rC(k)-rF(Nr+1) )/( rSigmaBnd-rF(Nr+1) )
281 c aHybSigmC(k) = bHybSigmC(k)*(rSigmaBnd-rF(1))*recip_fullDepth
282 c ENDIF
283 C- unique formula using max fct:
284 tmpRS = MAX( rC(k), rSigBndRS )
285 bHybSigmC(k) = ( rC(k)-rF(Nr+1) )/( tmpRS-rF(Nr+1) )
286 aHybSigmC(k) = bHybSigmC(k)*( tmpRS-rF(1) )*recip_fullDepth
287 ENDDO
288 ENDIF
289 ELSE
290 C-- Coeff at interface are read from file
291 IF (setCenterDr) THEN
292 STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID: Missing Code'
293 ENDIF
294 ENDIF
295
296 C-- Finish setting (if not done) using simple averaging
297 IF ( .NOT.setInterFDr ) THEN
298 C- Interface position at the middle between 2 centers
299 bHybSigmF(1) = 1. _d 0
300 aHybSigmF(1) = 0. _d 0
301 bHybSigmF(Nr+1) = 0. _d 0
302 aHybSigmF(Nr+1) = 0. _d 0
303 DO k=2,Nr
304 bHybSigmF(k) = ( bHybSigmC(k) + bHybSigmC(k-1) )*0.5 _d 0
305 aHybSigmF(k) = ( aHybSigmC(k) + aHybSigmC(k-1) )*0.5 _d 0
306 ENDDO
307 ENDIF
308 IF ( .NOT.setCenterDr ) THEN
309 C- Center position at the middle between 2 interfaces
310 DO k=1,Nr
311 bHybSigmC(k) = ( bHybSigmF(k) + bHybSigmF(k+1) )*0.5 _d 0
312 aHybSigmC(k) = ( aHybSigmF(k) + aHybSigmF(k+1) )*0.5 _d 0
313 ENDDO
314 ENDIF
315
316 C-- Vertical increment:
317 DO k=1,Nr
318 dAHybSigF(k) = ( aHybSigmF(k+1) - aHybSigmF(k) )*rkSign
319 dBHybSigF(k) = ( bHybSigmF(k+1) - bHybSigmF(k) )*rkSign
320 ENDDO
321 DO k=2,Nr
322 dAHybSigC(k) = ( aHybSigmC(k) - aHybSigmC(k-1) )*rkSign
323 dBHybSigC(k) = ( bHybSigmC(k) - bHybSigmC(k-1) )*rkSign
324 ENDDO
325 dAHybSigC(1) = ( aHybSigmC(1) - aHybSigmF(1) )*rkSign
326 dBHybSigC(1) = ( bHybSigmC(1) - bHybSigmF(1) )*rkSign
327 dAHybSigC(Nr+1) = ( aHybSigmF(Nr+1) - aHybSigmC(Nr) )*rkSign
328 dBHybSigC(Nr+1) = ( bHybSigmF(Nr+1) - bHybSigmC(Nr) )*rkSign
329
330 C-- Check for miss-match between vertical discretization :
331 maxErrC = 0.
332 maxErrF = 0.
333 epsil = 1. _d -9
334 DO k=1,Nr
335 tmpError = ( rC(k)-rF(Nr+1) )*recip_fullDepth
336 & - ( aHybSigmC(k)+bHybSigmC(k) )
337 IF ( ABS(tmpError).GT.epsil ) THEN
338 IF ( maxErrC.LE.epsil ) THEN
339 WRITE(msgBuf,'(2A)') 'S/R INI_VERTICAL_GRID:',
340 & ' rC and Hybrid-Sigma Coeff miss-match'
341 CALL PRINT_ERROR( msgBuf, myThid )
342 ENDIF
343 WRITE(msgBuf,'(A,I4,2(A,1PE14.6),A,1P2E14.6)')
344 & ' k=', k,' , err=', tmpError, ' ; rC=', rC(k),
345 & ' , a & b=', aHybSigmC(k), bHybSigmC(k)
346 CALL PRINT_ERROR( msgBuf, myThid )
347 ENDIF
348 maxErrC = MAX( maxErrC, ABS(tmpError) )
349 ENDDO
350 DO k=1,Nr+1
351 tmpError = ( rF(k)-rF(Nr+1) )*recip_fullDepth
352 & - ( aHybSigmF(k)+bHybSigmF(k) )
353 IF ( ABS(tmpError).GT.epsil ) THEN
354 IF ( maxErrF.LE.epsil ) THEN
355 WRITE(msgBuf,'(2A)') 'S/R INI_VERTICAL_GRID:',
356 & ' rF and Hybrid-Sigma Coeff miss-match'
357 CALL PRINT_ERROR( msgBuf, myThid )
358 ENDIF
359 WRITE(msgBuf,'(A,I4,2(A,1PE14.6),A,1P2E14.6)')
360 & ' k=', k,' , err=', tmpError, ' ; rF=', rF(k),
361 & ' , a & b=', aHybSigmF(k), bHybSigmF(k)
362 CALL PRINT_ERROR( msgBuf, myThid )
363 ENDIF
364 maxErrF = MAX( maxErrF, ABS(tmpError) )
365 ENDDO
366 WRITE(msgBuf,'(2A,1PE14.6)') 'S/R INI_VERTICAL_GRID:',
367 & ' matching of aHybSigmC & rC :', maxErrC
368 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
369 & SQUEEZE_RIGHT, myThid )
370 WRITE(msgBuf,'(2A,1PE14.6)') 'S/R INI_VERTICAL_GRID:',
371 & ' matching of aHybSigmF & rF :', maxErrF
372 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
373 & SQUEEZE_RIGHT, myThid )
374 IF ( maxErrC.GT.epsil .OR. maxErrF.GT.epsil ) THEN
375 WRITE(msgBuf,'(2A)') 'S/R INI_VERTICAL_GRID:',
376 & ' rC,rF and Hybrid-Sigma Coeff miss-match'
377 CALL PRINT_ERROR( msgBuf, myThid )
378 STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID'
379 ENDIF
380
381 C--- End setting-up Hybrid-Sigma vertical coordinate:
382 ENDIF
383
384 _END_MASTER(myThid)
385 _BARRIER
386
387 RETURN
388 END

  ViewVC Help
Powered by ViewVC 1.1.22