/[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.23 - (show annotations) (download)
Mon Apr 4 21:29:00 2016 UTC (8 years ago) by jmc
Branch: MAIN
CVS Tags: checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65v, checkpoint65w, checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, HEAD
Changes since 1.22: +39 -9 lines
- new parameters "top_Pres" & "seaLev_Z" (replacing Ro_SeaLevel and recently
  added phi0Ref) to set vertical axis origin and phiRef origin according to
  coordinate and fluid type.

1 C $Header: /u/gcmpack/MITgcm/model/src/ini_vertical_grid.F,v 1.22 2013/07/26 14:22:22 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 C -- Use top_Pres or seaLev_Z to set vertical axis position:
135 C OCN in Z: top_Pres(Ref) (= rhoConst*PhiRef(1)), seaLev_Z (=rF(1), @ the top)
136 C ATM in Z: top_Pres(Ref) (= rhoConst*PhiRef(1)), seaLev_Z (=rF(Nr+1), @ bottom)
137 C OCN in P: top_Pres (=rF(Nr+1)), seaLev_Z (<-> PhiRef(Nr+1)/g, @ the top)
138 C ATM in P: top_Pres (=rF(Nr+1)), seaLev_Z (<-> PhiRef(1)/g, @ the bottom)
139 IF ( rF(1).EQ.UNSET_RS .AND.
140 & usingZCoords.AND.fluidIsWater ) THEN
141 rF(1) = seaLev_Z
142 ENDIF
143 IF ( rF(1).NE.UNSET_RS ) THEN
144 DO k=1,Nr
145 rF(k+1) = rF(k) + rkSign*drF(k)
146 ENDDO
147 rC(1) = rF(1) + rkSign*drC(1)
148 DO k=2,Nr
149 rC(k) = rC(k-1) + rkSign*drC(k)
150 ENDDO
151 c IF ( usingPCoords ) THEN
152 c top_Pres = rF(Nr+1)
153 c ELSEIF ( fluidIsWater ) THEN
154 c seaLev_Z = rF(1)
155 c ELSE
156 c seaLev_Z = rF(Nr+1)
157 c ENDIF
158 ELSE
159 IF ( usingPCoords ) THEN
160 rF(Nr+1) = top_Pres
161 ELSE
162 rF(Nr+1) = seaLev_Z
163 ENDIF
164 DO k=Nr,1,-1
165 rF(k) = rF(k+1) - rkSign*drF(k)
166 ENDDO
167 rC(Nr) = rF(Nr+1) - rkSign*drC(Nr+1)
168 DO k=Nr,2,-1
169 rC(k-1) = rC(k) - rkSign*drC(k)
170 ENDDO
171 ENDIF
172
173 C--- Check vertical discretization :
174 checkRatio2 = 100.
175 checkRatio1 = 1. _d 0 / checkRatio2
176 DO k=1,Nr
177 tmpRatio = 0.
178 IF ( (rC(k)-rF(k+1)) .NE. 0. )
179 & tmpRatio = (rF(k)-rC(k)) / (rC(k)-rF(k+1))
180 IF ( tmpRatio.LT.checkRatio1 .OR. tmpRatio.GT.checkRatio2 ) THEN
181 c write(0,*) 'drF=',drF
182 c write(0,*) 'drC=',drC
183 c write(0,*) 'rF=',rF
184 c write(0,*) 'rC=',rC
185 WRITE(msgBuf,'(A,I4,A,E16.8)')
186 & 'S/R INI_VERTICAL_GRID: Invalid relative position, level k=',
187 & k, ' :', tmpRatio
188 CALL PRINT_ERROR( msgBuf, myThid )
189 WRITE(msgBuf,'(A,1PE14.6,A,2E14.6)')
190 & 'S/R INI_VERTICAL_GRID: rC=', rC(k),
191 & ' , rF(k,k+1)=',rF(k),rF(k+1)
192 CALL PRINT_ERROR( msgBuf, myThid )
193 STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID'
194 ENDIF
195 ENDDO
196
197 C- Calculate reciprol vertical grid spacing :
198 DO k=1,Nr+1
199 recip_drC(k) = 1. _d 0/drC(k)
200 ENDDO
201 DO k=1,Nr
202 recip_drF(k) = 1. _d 0/drF(k)
203 ENDDO
204
205 C--- Hybrid-Sigma vertical coordinate:
206 IF ( selectSigmaCoord .EQ. 0 ) THEN
207 DO k=1,Nr+1
208 aHybSigmF(k) = 0. _d 0
209 bHybSigmF(k) = 0. _d 0
210 dAHybSigC(k) = 0. _d 0
211 dAHybSigC(k) = 0. _d 0
212 ENDDO
213 DO k=1,Nr
214 aHybSigmC(k) = 0. _d 0
215 bHybSigmC(k) = 0. _d 0
216 dAHybSigF(k) = 0. _d 0
217 dAHybSigF(k) = 0. _d 0
218 ENDDO
219 ELSE
220 rFullDepth = rF(1) - rF(Nr+1)
221 recip_fullDepth = 0. _d 0
222 IF ( rFullDepth.GT.0. ) recip_fullDepth = 1. _d 0 / rFullDepth
223 rSigBndRS = rSigmaBnd
224 IF ( hybSigmFile.EQ.' ' .AND. rSigmaBnd.EQ.UNSET_RL ) THEN
225 C- Default is pure sigma:
226 IF ( usingPCoords ) rSigBndRS = rF(Nr+1)
227 IF ( usingZCoords ) rSigBndRS = rF(1)
228 ENDIF
229 c IF ( hybSigmFile.EQ.' ' .AND. rSigmaBnd.EQ.UNSET_RL ) THEN
230 C- compute coeff as pure sigma coordinate
231 c DO k=1,Nr+1
232 c aHybSigmF(k) = 0.
233 c bHybSigmF(k) = (rF(k)-rF(Nr+1))*recip_fullDepth
234 c ENDDO
235 c DO k=1,Nr
236 c aHybSigmC(k) = 0.
237 c bHybSigmC(k) = (rC(k)-rF(Nr+1))*recip_fullDepth
238 c ENDDO
239 c ELSEIF ( hybSigmFile.EQ.' ' ) THEN
240 IF ( hybSigmFile.EQ.' ' ) THEN
241 C-- compute coeff assuming fixed r-coord above rSigmaBnd and pure sigma below
242 IF ( usingPCoords .AND. setInterFDr ) THEN
243 C- Alternative method : p = pTop + A*DeltaFullP + B*(eta+Pr_surf - rSigmaBnd)
244 c aHybSigmF(k) = ( MIN(rF(k),rSigmaBnd) - rF(Nr+1) )
245 c & *recip_fullDepth
246 c bHybSigmF(k) = ( MAX(rF(k),rSigmaBnd) - rSigmaBnd )
247 c & /( rF(1) - rSigmaBnd )
248 C- Standard method : p = pTop + A*DeltaFullP + B*(eta+Ro_surf - pTop)
249 C sigma part goes from 0 @ rSigmaBnd (and above) to 1 @ surface
250 DO k=1,Nr+1
251 C- separate the 2 cases:
252 c IF ( rF(k).LE.rSigmaBnd ) THEN
253 c bHybSigmF(k) = 0.
254 c aHybSigmF(k) = ( rF(k) - rF(Nr+1) )*recip_fullDepth
255 c ELSE
256 c bHybSigmF(k) = (rF(k)-rSigmaBnd)/(rF(1)-rSigmaBnd)
257 c aHybSigmF(k) = (1. _d 0 - bHybSigmF(k))
258 c & *(rSigmaBnd-rF(Nr+1) )*recip_fullDepth
259 c ENDIF
260 C- unique formula using min fct:
261 tmpRS = MIN( rF(k), rSigBndRS )
262 bHybSigmF(k) = ( rF(k) - tmpRS )/(rF(1)-rSigBndRS)
263 aHybSigmF(k) = (1. _d 0 - bHybSigmF(k))
264 & *( tmpRS -rF(Nr+1) )*recip_fullDepth
265 ENDDO
266 ENDIF
267 IF ( usingPCoords .AND. setCenterDr ) THEN
268 DO k=1,Nr
269 C- separate the 2 cases:
270 c IF ( rC(k).LE.rSigmaBnd ) THEN
271 c bHybSigmC(k) = 0.
272 c aHybSigmC(k) = ( rC(k) - rF(Nr+1) )*recip_fullDepth
273 c ELSE
274 c bHybSigmC(k) = (rC(k)-rSigmaBnd)/(rF(1)-rSigmaBnd)
275 c aHybSigmC(k) = (1. _d 0 - bHybSigmC(k))
276 c & *(rSigmaBnd-rF(Nr+1) )*recip_fullDepth
277 c ENDIF
278 C- unique formula using min fct:
279 tmpRS = MIN( rC(k), rSigBndRS )
280 bHybSigmC(k) = ( rC(k) - tmpRS )/(rF(1)-rSigBndRS)
281 aHybSigmC(k) = (1. _d 0 - bHybSigmC(k))
282 & *( tmpRS -rF(Nr+1) )*recip_fullDepth
283 ENDDO
284 ENDIF
285 IF ( usingZCoords .AND. setInterFDr ) THEN
286 C- Standard method : z = zBot + A*DeltaFullZ + B*(eta+Ro_surf - zBot)
287 C sigma part goes from 1 @ rSigmaBnd (and above) to 0 @ bottom
288 DO k=1,Nr+1
289 C- separate the 2 cases:
290 c IF ( rF(k).GE.rSigmaBnd ) THEN
291 c bHybSigmF(k) = 1.
292 c aHybSigmF(k) = ( rF(k)-rF(1) )*recip_fullDepth
293 c ELSE
294 c bHybSigmF(k) = ( rF(k)-rF(Nr+1) )/( rSigmaBnd-rF(Nr+1) )
295 c aHybSigmF(k) = bHybSigmF(k)*(rSigmaBnd-rF(1))*recip_fullDepth
296 c ENDIF
297 C- unique formula using max fct:
298 tmpRS = MAX( rF(k), rSigBndRS )
299 bHybSigmF(k) = ( rF(k)-rF(Nr+1) )/( tmpRS-rF(Nr+1) )
300 aHybSigmF(k) = bHybSigmF(k)*( tmpRS-rF(1) )*recip_fullDepth
301 ENDDO
302 ENDIF
303 IF ( usingZCoords .AND. setCenterDr ) THEN
304 DO k=1,Nr
305 C- separate the 2 cases:
306 c IF ( rC(k).GE.rSigmaBnd ) THEN
307 c bHybSigmC(k) = 1.
308 c aHybSigmC(k) = ( rC(k)-rF(1) )*recip_fullDepth
309 c ELSE
310 c bHybSigmC(k) = ( rC(k)-rF(Nr+1) )/( rSigmaBnd-rF(Nr+1) )
311 c aHybSigmC(k) = bHybSigmC(k)*(rSigmaBnd-rF(1))*recip_fullDepth
312 c ENDIF
313 C- unique formula using max fct:
314 tmpRS = MAX( rC(k), rSigBndRS )
315 bHybSigmC(k) = ( rC(k)-rF(Nr+1) )/( tmpRS-rF(Nr+1) )
316 aHybSigmC(k) = bHybSigmC(k)*( tmpRS-rF(1) )*recip_fullDepth
317 ENDDO
318 ENDIF
319 ELSE
320 C-- Coeff at interface are read from file
321 IF (setCenterDr) THEN
322 STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID: Missing Code'
323 ENDIF
324 ENDIF
325
326 C-- Finish setting (if not done) using simple averaging
327 IF ( .NOT.setInterFDr ) THEN
328 C- Interface position at the middle between 2 centers
329 bHybSigmF(1) = 1. _d 0
330 aHybSigmF(1) = 0. _d 0
331 bHybSigmF(Nr+1) = 0. _d 0
332 aHybSigmF(Nr+1) = 0. _d 0
333 DO k=2,Nr
334 bHybSigmF(k) = ( bHybSigmC(k) + bHybSigmC(k-1) )*0.5 _d 0
335 aHybSigmF(k) = ( aHybSigmC(k) + aHybSigmC(k-1) )*0.5 _d 0
336 ENDDO
337 ENDIF
338 IF ( .NOT.setCenterDr ) THEN
339 C- Center position at the middle between 2 interfaces
340 DO k=1,Nr
341 bHybSigmC(k) = ( bHybSigmF(k) + bHybSigmF(k+1) )*0.5 _d 0
342 aHybSigmC(k) = ( aHybSigmF(k) + aHybSigmF(k+1) )*0.5 _d 0
343 ENDDO
344 ENDIF
345
346 C-- Vertical increment:
347 DO k=1,Nr
348 dAHybSigF(k) = ( aHybSigmF(k+1) - aHybSigmF(k) )*rkSign
349 dBHybSigF(k) = ( bHybSigmF(k+1) - bHybSigmF(k) )*rkSign
350 ENDDO
351 DO k=2,Nr
352 dAHybSigC(k) = ( aHybSigmC(k) - aHybSigmC(k-1) )*rkSign
353 dBHybSigC(k) = ( bHybSigmC(k) - bHybSigmC(k-1) )*rkSign
354 ENDDO
355 dAHybSigC(1) = ( aHybSigmC(1) - aHybSigmF(1) )*rkSign
356 dBHybSigC(1) = ( bHybSigmC(1) - bHybSigmF(1) )*rkSign
357 dAHybSigC(Nr+1) = ( aHybSigmF(Nr+1) - aHybSigmC(Nr) )*rkSign
358 dBHybSigC(Nr+1) = ( bHybSigmF(Nr+1) - bHybSigmC(Nr) )*rkSign
359
360 C-- Check for miss-match between vertical discretization :
361 maxErrC = 0.
362 maxErrF = 0.
363 epsil = 1. _d -9
364 DO k=1,Nr
365 tmpError = ( rC(k)-rF(Nr+1) )*recip_fullDepth
366 & - ( aHybSigmC(k)+bHybSigmC(k) )
367 IF ( ABS(tmpError).GT.epsil ) THEN
368 IF ( maxErrC.LE.epsil ) THEN
369 WRITE(msgBuf,'(2A)') 'S/R INI_VERTICAL_GRID:',
370 & ' rC and Hybrid-Sigma Coeff miss-match'
371 CALL PRINT_ERROR( msgBuf, myThid )
372 ENDIF
373 WRITE(msgBuf,'(A,I4,2(A,1PE14.6),A,1P2E14.6)')
374 & ' k=', k,' , err=', tmpError, ' ; rC=', rC(k),
375 & ' , a & b=', aHybSigmC(k), bHybSigmC(k)
376 CALL PRINT_ERROR( msgBuf, myThid )
377 ENDIF
378 maxErrC = MAX( maxErrC, ABS(tmpError) )
379 ENDDO
380 DO k=1,Nr+1
381 tmpError = ( rF(k)-rF(Nr+1) )*recip_fullDepth
382 & - ( aHybSigmF(k)+bHybSigmF(k) )
383 IF ( ABS(tmpError).GT.epsil ) THEN
384 IF ( maxErrF.LE.epsil ) THEN
385 WRITE(msgBuf,'(2A)') 'S/R INI_VERTICAL_GRID:',
386 & ' rF and Hybrid-Sigma Coeff miss-match'
387 CALL PRINT_ERROR( msgBuf, myThid )
388 ENDIF
389 WRITE(msgBuf,'(A,I4,2(A,1PE14.6),A,1P2E14.6)')
390 & ' k=', k,' , err=', tmpError, ' ; rF=', rF(k),
391 & ' , a & b=', aHybSigmF(k), bHybSigmF(k)
392 CALL PRINT_ERROR( msgBuf, myThid )
393 ENDIF
394 maxErrF = MAX( maxErrF, ABS(tmpError) )
395 ENDDO
396 WRITE(msgBuf,'(2A,1PE14.6)') 'S/R INI_VERTICAL_GRID:',
397 & ' matching of aHybSigmC & rC :', maxErrC
398 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
399 & SQUEEZE_RIGHT, myThid )
400 WRITE(msgBuf,'(2A,1PE14.6)') 'S/R INI_VERTICAL_GRID:',
401 & ' matching of aHybSigmF & rF :', maxErrF
402 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
403 & SQUEEZE_RIGHT, myThid )
404 IF ( maxErrC.GT.epsil .OR. maxErrF.GT.epsil ) THEN
405 WRITE(msgBuf,'(2A)') 'S/R INI_VERTICAL_GRID:',
406 & ' rC,rF and Hybrid-Sigma Coeff miss-match'
407 CALL PRINT_ERROR( msgBuf, myThid )
408 STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID'
409 ENDIF
410
411 C--- End setting-up Hybrid-Sigma vertical coordinate:
412 ENDIF
413
414 _END_MASTER(myThid)
415 _BARRIER
416
417 RETURN
418 END

  ViewVC Help
Powered by ViewVC 1.1.22