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

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

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


Revision 1.22 - (hide annotations) (download)
Fri Jul 26 14:22:22 2013 UTC (10 years, 9 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 jmc 1.22 C $Header: /u/gcmpack/MITgcm/model/src/ini_vertical_grid.F,v 1.21 2011/08/18 20:29:49 jmc Exp $
2 cnh 1.12 C $Name: $
3 adcroft 1.1
4 cnh 1.7 #include "CPP_OPTIONS.h"
5 adcroft 1.1
6 cnh 1.12 CBOP
7     C !ROUTINE: INI_VERTICAL_GRID
8     C !INTERFACE:
9 adcroft 1.1 SUBROUTINE INI_VERTICAL_GRID( myThid )
10 jmc 1.17
11 cnh 1.12 C !DESCRIPTION: \bv
12     C *==========================================================*
13 jmc 1.17 C | SUBROUTINE INI_VERTICAL_GRID
14     C | o Initialise vertical gridding arrays
15 cnh 1.12 C *==========================================================*
16     C \ev
17    
18     C !USES:
19 adcroft 1.8 IMPLICIT NONE
20 adcroft 1.1 C === Global variables ===
21     #include "SIZE.h"
22     #include "EEPARAMS.h"
23     #include "PARAMS.h"
24     #include "GRID.h"
25    
26 cnh 1.12 C !INPUT/OUTPUT PARAMETERS:
27 adcroft 1.1 C == Routine arguments ==
28 jmc 1.17 C myThid :: my Thread Id number
29 adcroft 1.1 INTEGER myThid
30    
31 cnh 1.12 C !LOCAL VARIABLES:
32 adcroft 1.1 C == Local variables ==
33 jmc 1.17 C k :: loop index
34 jmc 1.20 C msgBuf :: Informational/error message buffer
35 jmc 1.17 INTEGER k
36     _RL tmpRatio, checkRatio1, checkRatio2
37 jmc 1.13 CHARACTER*(MAX_LEN_MBUF) msgBuf
38 jmc 1.20 _RL maxErrC, maxErrF, epsil, tmpError
39     _RL rFullDepth, recip_fullDepth
40     _RS rSigBndRS, tmpRS
41 cnh 1.12 CEOP
42 adcroft 1.1
43 jmc 1.15 _BEGIN_MASTER(myThid)
44    
45 jmc 1.18 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 jmc 1.16 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 jmc 1.17 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 jmc 1.16 rkSign = -1. _d 0
57     gravitySign = -1. _d 0
58     IF ( usingPCoords ) THEN
59     gravitySign = 1. _d 0
60     ENDIF
61    
62 jmc 1.18 IF ( .NOT.(setInterFDr.OR.setCenterDr) ) THEN
63     WRITE(msgBuf,'(A)')
64 jmc 1.17 & 'S/R INI_VERTICAL_GRID: neither delR nor delRc are defined'
65 jmc 1.18 CALL PRINT_ERROR( msgBuf, myThid )
66     WRITE(msgBuf,'(A)')
67 jmc 1.17 & 'S/R INI_VERTICAL_GRID: Need at least 1 of the 2 (delR,delRc)'
68 jmc 1.18 CALL PRINT_ERROR( msgBuf, myThid )
69     STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID'
70 jmc 1.17 ENDIF
71 jmc 1.13
72 jmc 1.17 C--- Set Level r-thickness (drF) and Center r-distances (drC)
73 jmc 1.18
74 jmc 1.17 IF (setInterFDr) THEN
75     C-- Interface r-distances are defined:
76     DO k=1,Nr
77 jmc 1.18 drF(k) = delR(k)
78 jmc 1.17 ENDDO
79 jmc 1.13 C- Check that all thickness are > 0 :
80 jmc 1.17 DO k=1,Nr
81 jmc 1.18 IF (delR(k).LE.0.) THEN
82     WRITE(msgBuf,'(A,I4,A,E16.8)')
83 jmc 1.17 & 'S/R INI_VERTICAL_GRID: delR(k=',k,' )=',delR(k)
84 jmc 1.18 CALL PRINT_ERROR( msgBuf, myThid )
85     WRITE(msgBuf,'(A)')
86 jmc 1.13 & 'S/R INI_VERTICAL_GRID: Vert. grid spacing MUST BE > 0'
87 jmc 1.18 CALL PRINT_ERROR( msgBuf, myThid )
88     STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID'
89     ENDIF
90 jmc 1.17 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 jmc 1.21 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 jmc 1.18 drF(k-1) = 0.5 _d 0 *delRc(k) + drF(k-1)
102 jmc 1.13 ENDDO
103 jmc 1.17 drF(Nr) = delRc(Nr+1) + drF(Nr)
104     ENDIF
105 jmc 1.13
106 jmc 1.17 IF (setCenterDr) THEN
107     C-- Cell Center r-distances are defined:
108 jmc 1.22 DO k=1,Nr+1
109 jmc 1.17 drC(k) = delRc(k)
110 jmc 1.13 ENDDO
111     C- Check that all thickness are > 0 :
112 jmc 1.17 DO k=1,Nr+1
113     IF (delRc(k).LE.0.) THEN
114 jmc 1.13 WRITE(msgBuf,'(A,I4,A,E16.8)')
115 jmc 1.17 & 'S/R INI_VERTICAL_GRID: delRc(k=',k,' )=',delRc(k)
116     CALL PRINT_ERROR( msgBuf, myThid )
117 jmc 1.13 WRITE(msgBuf,'(A)')
118     & 'S/R INI_VERTICAL_GRID: Vert. grid spacing MUST BE > 0'
119 jmc 1.17 CALL PRINT_ERROR( msgBuf, myThid )
120 jmc 1.13 STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID'
121     ENDIF
122     ENDDO
123 jmc 1.17 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 jmc 1.22 drC(Nr+1) = 0.5 _d 0 *delR(Nr)
131 jmc 1.17 ENDIF
132 jmc 1.13
133 jmc 1.17 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 adcroft 1.1 ENDDO
138 jmc 1.17 rC(1) = rF(1) + rkSign*drC(1)
139     DO k=2,Nr
140     rC(k) = rC(k-1) + rkSign*drC(k)
141 adcroft 1.1 ENDDO
142 jmc 1.13
143 jmc 1.17 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 jmc 1.13
167     C- Calculate reciprol vertical grid spacing :
168 jmc 1.22 DO k=1,Nr+1
169     recip_drC(k) = 1. _d 0/drC(k)
170     ENDDO
171 jmc 1.17 DO k=1,Nr
172     recip_drF(k) = 1. _d 0/drF(k)
173 adcroft 1.1 ENDDO
174 jmc 1.13
175 jmc 1.20 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 jmc 1.15 _END_MASTER(myThid)
385     _BARRIER
386    
387 adcroft 1.1 RETURN
388     END

  ViewVC Help
Powered by ViewVC 1.1.22