/[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.20 - (hide annotations) (download)
Sat Sep 11 21:24:52 2010 UTC (13 years, 8 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint62v, checkpoint62u, checkpoint62t, checkpoint62s, checkpoint62r, checkpoint62q, checkpoint62p, checkpoint62k, checkpoint62o, checkpoint62n, checkpoint62m, checkpoint62l, checkpoint62w, checkpoint62z, checkpoint62y, checkpoint62x, checkpoint63, checkpoint63a
Changes since 1.19: +214 -2 lines
first check-in of sigma (and hybrid-sigma) coordinate code

1 jmc 1.20 C $Header: /u/gcmpack/MITgcm/model/src/ini_vertical_grid.F,v 1.19 2008/09/05 20:15:28 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.18 drF(k-1) = 0.5 _d 0 *delRc(k) + drF(k-1)
97     drF( k ) = 0.5 _d 0 *delRc(k)
98 jmc 1.13 ENDDO
99 jmc 1.17 drF(Nr) = delRc(Nr+1) + drF(Nr)
100     ENDIF
101 jmc 1.13
102 jmc 1.17 IF (setCenterDr) THEN
103     C-- Cell Center r-distances are defined:
104     DO k=1,Nr
105     drC(k) = delRc(k)
106 jmc 1.13 ENDDO
107     C- Check that all thickness are > 0 :
108 jmc 1.17 DO k=1,Nr+1
109     IF (delRc(k).LE.0.) THEN
110 jmc 1.13 WRITE(msgBuf,'(A,I4,A,E16.8)')
111 jmc 1.17 & 'S/R INI_VERTICAL_GRID: delRc(k=',k,' )=',delRc(k)
112     CALL PRINT_ERROR( msgBuf, myThid )
113 jmc 1.13 WRITE(msgBuf,'(A)')
114     & 'S/R INI_VERTICAL_GRID: Vert. grid spacing MUST BE > 0'
115 jmc 1.17 CALL PRINT_ERROR( msgBuf, myThid )
116 jmc 1.13 STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID'
117     ENDIF
118     ENDDO
119 jmc 1.17 ELSE
120     C-- Cell Center r-distances undefined:
121     C assume Center at middle between 2 Interfaces
122     drC(1) = 0.5 _d 0 *delR(1)
123     DO k=2,Nr
124     drC(k) = 0.5 _d 0 *(delR(k-1)+delR(k))
125     ENDDO
126     ENDIF
127 jmc 1.13
128 jmc 1.17 C--- Set r-position of interFace (rF) and cell-Center (rC):
129     rF(1) = Ro_SeaLevel
130     DO k=1,Nr
131     rF(k+1) = rF(k) + rkSign*drF(k)
132 adcroft 1.1 ENDDO
133 jmc 1.17 rC(1) = rF(1) + rkSign*drC(1)
134     DO k=2,Nr
135     rC(k) = rC(k-1) + rkSign*drC(k)
136 adcroft 1.1 ENDDO
137 jmc 1.13
138 jmc 1.17 C--- Check vertical discretization :
139     checkRatio2 = 100.
140     checkRatio1 = 1. _d 0 / checkRatio2
141     DO k=1,Nr
142     tmpRatio = 0.
143     IF ( (rC(k)-rF(k+1)) .NE. 0. )
144     & tmpRatio = (rF(k)-rC(k)) / (rC(k)-rF(k+1))
145     IF ( tmpRatio.LT.checkRatio1 .OR. tmpRatio.GT.checkRatio2 ) THEN
146     c write(0,*) 'drF=',drF
147     c write(0,*) 'drC=',drC
148     c write(0,*) 'rF=',rF
149     c write(0,*) 'rC=',rC
150     WRITE(msgBuf,'(A,I4,A,E16.8)')
151     & 'S/R INI_VERTICAL_GRID: Invalid relative position, level k=',
152     & k, ' :', tmpRatio
153     CALL PRINT_ERROR( msgBuf, myThid )
154     WRITE(msgBuf,'(A,1PE14.6,A,2E14.6)')
155     & 'S/R INI_VERTICAL_GRID: rC=', rC(k),
156     & ' , rF(k,k+1)=',rF(k),rF(k+1)
157     CALL PRINT_ERROR( msgBuf, myThid )
158     STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID'
159     ENDIF
160     ENDDO
161 jmc 1.13
162     C- Calculate reciprol vertical grid spacing :
163 jmc 1.17 DO k=1,Nr
164     recip_drC(k) = 1. _d 0/drC(k)
165     recip_drF(k) = 1. _d 0/drF(k)
166 adcroft 1.1 ENDDO
167 jmc 1.13
168 jmc 1.20 C--- Hybrid-Sigma vertical coordinate:
169     IF ( selectSigmaCoord .EQ. 0 ) THEN
170     DO k=1,Nr+1
171     aHybSigmF(k) = 0. _d 0
172     bHybSigmF(k) = 0. _d 0
173     dAHybSigC(k) = 0. _d 0
174     dAHybSigC(k) = 0. _d 0
175     ENDDO
176     DO k=1,Nr
177     aHybSigmC(k) = 0. _d 0
178     bHybSigmC(k) = 0. _d 0
179     dAHybSigF(k) = 0. _d 0
180     dAHybSigF(k) = 0. _d 0
181     ENDDO
182     ELSE
183     rFullDepth = rF(1) - rF(Nr+1)
184     recip_fullDepth = 0. _d 0
185     IF ( rFullDepth.GT.0. ) recip_fullDepth = 1. _d 0 / rFullDepth
186     rSigBndRS = rSigmaBnd
187     IF ( hybSigmFile.EQ.' ' .AND. rSigmaBnd.EQ.UNSET_RL ) THEN
188     C- Default is pure sigma:
189     IF ( usingPCoords ) rSigBndRS = rF(Nr+1)
190     IF ( usingZCoords ) rSigBndRS = rF(1)
191     ENDIF
192     c IF ( hybSigmFile.EQ.' ' .AND. rSigmaBnd.EQ.UNSET_RL ) THEN
193     C- compute coeff as pure sigma coordinate
194     c DO k=1,Nr+1
195     c aHybSigmF(k) = 0.
196     c bHybSigmF(k) = (rF(k)-rF(Nr+1))*recip_fullDepth
197     c ENDDO
198     c DO k=1,Nr
199     c aHybSigmC(k) = 0.
200     c bHybSigmC(k) = (rC(k)-rF(Nr+1))*recip_fullDepth
201     c ENDDO
202     c ELSEIF ( hybSigmFile.EQ.' ' ) THEN
203     IF ( hybSigmFile.EQ.' ' ) THEN
204     C-- compute coeff assuming fixed r-coord above rSigmaBnd and pure sigma below
205     IF ( usingPCoords .AND. setInterFDr ) THEN
206     C- Alternative method : p = pTop + A*DeltaFullP + B*(eta+Pr_surf - rSigmaBnd)
207     c aHybSigmF(k) = ( MIN(rF(k),rSigmaBnd) - rF(Nr+1) )
208     c & *recip_fullDepth
209     c bHybSigmF(k) = ( MAX(rF(k),rSigmaBnd) - rSigmaBnd )
210     c & /( rF(1) - rSigmaBnd )
211     C- Standard method : p = pTop + A*DeltaFullP + B*(eta+Ro_surf - pTop)
212     C sigma part goes from 0 @ rSigmaBnd (and above) to 1 @ surface
213     DO k=1,Nr+1
214     C- separate the 2 cases:
215     c IF ( rF(k).LE.rSigmaBnd ) THEN
216     c bHybSigmF(k) = 0.
217     c aHybSigmF(k) = ( rF(k) - rF(Nr+1) )*recip_fullDepth
218     c ELSE
219     c bHybSigmF(k) = (rF(k)-rSigmaBnd)/(rF(1)-rSigmaBnd)
220     c aHybSigmF(k) = (1. _d 0 - bHybSigmF(k))
221     c & *(rSigmaBnd-rF(Nr+1) )*recip_fullDepth
222     c ENDIF
223     C- unique formula using min fct:
224     tmpRS = MIN( rF(k), rSigBndRS )
225     bHybSigmF(k) = ( rF(k) - tmpRS )/(rF(1)-rSigBndRS)
226     aHybSigmF(k) = (1. _d 0 - bHybSigmF(k))
227     & *( tmpRS -rF(Nr+1) )*recip_fullDepth
228     ENDDO
229     ENDIF
230     IF ( usingPCoords .AND. setCenterDr ) THEN
231     DO k=1,Nr
232     C- separate the 2 cases:
233     c IF ( rC(k).LE.rSigmaBnd ) THEN
234     c bHybSigmC(k) = 0.
235     c aHybSigmC(k) = ( rC(k) - rF(Nr+1) )*recip_fullDepth
236     c ELSE
237     c bHybSigmC(k) = (rC(k)-rSigmaBnd)/(rF(1)-rSigmaBnd)
238     c aHybSigmC(k) = (1. _d 0 - bHybSigmC(k))
239     c & *(rSigmaBnd-rF(Nr+1) )*recip_fullDepth
240     c ENDIF
241     C- unique formula using min fct:
242     tmpRS = MIN( rC(k), rSigBndRS )
243     bHybSigmC(k) = ( rC(k) - tmpRS )/(rF(1)-rSigBndRS)
244     aHybSigmC(k) = (1. _d 0 - bHybSigmC(k))
245     & *( tmpRS -rF(Nr+1) )*recip_fullDepth
246     ENDDO
247     ENDIF
248     IF ( usingZCoords .AND. setInterFDr ) THEN
249     C- Standard method : z = zBot + A*DeltaFullZ + B*(eta+Ro_surf - zBot)
250     C sigma part goes from 1 @ rSigmaBnd (and above) to 0 @ bottom
251     DO k=1,Nr+1
252     C- separate the 2 cases:
253     c IF ( rF(k).GE.rSigmaBnd ) THEN
254     c bHybSigmF(k) = 1.
255     c aHybSigmF(k) = ( rF(k)-rF(1) )*recip_fullDepth
256     c ELSE
257     c bHybSigmF(k) = ( rF(k)-rF(Nr+1) )/( rSigmaBnd-rF(Nr+1) )
258     c aHybSigmF(k) = bHybSigmF(k)*(rSigmaBnd-rF(1))*recip_fullDepth
259     c ENDIF
260     C- unique formula using max fct:
261     tmpRS = MAX( rF(k), rSigBndRS )
262     bHybSigmF(k) = ( rF(k)-rF(Nr+1) )/( tmpRS-rF(Nr+1) )
263     aHybSigmF(k) = bHybSigmF(k)*( tmpRS-rF(1) )*recip_fullDepth
264     ENDDO
265     ENDIF
266     IF ( usingZCoords .AND. setCenterDr ) THEN
267     DO k=1,Nr
268     C- separate the 2 cases:
269     c IF ( rC(k).GE.rSigmaBnd ) THEN
270     c bHybSigmC(k) = 1.
271     c aHybSigmC(k) = ( rC(k)-rF(1) )*recip_fullDepth
272     c ELSE
273     c bHybSigmC(k) = ( rC(k)-rF(Nr+1) )/( rSigmaBnd-rF(Nr+1) )
274     c aHybSigmC(k) = bHybSigmC(k)*(rSigmaBnd-rF(1))*recip_fullDepth
275     c ENDIF
276     C- unique formula using max fct:
277     tmpRS = MAX( rC(k), rSigBndRS )
278     bHybSigmC(k) = ( rC(k)-rF(Nr+1) )/( tmpRS-rF(Nr+1) )
279     aHybSigmC(k) = bHybSigmC(k)*( tmpRS-rF(1) )*recip_fullDepth
280     ENDDO
281     ENDIF
282     ELSE
283     C-- Coeff at interface are read from file
284     IF (setCenterDr) THEN
285     STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID: Missing Code'
286     ENDIF
287     ENDIF
288    
289     C-- Finish setting (if not done) using simple averaging
290     IF ( .NOT.setInterFDr ) THEN
291     C- Interface position at the middle between 2 centers
292     bHybSigmF(1) = 1. _d 0
293     aHybSigmF(1) = 0. _d 0
294     bHybSigmF(Nr+1) = 0. _d 0
295     aHybSigmF(Nr+1) = 0. _d 0
296     DO k=2,Nr
297     bHybSigmF(k) = ( bHybSigmC(k) + bHybSigmC(k-1) )*0.5 _d 0
298     aHybSigmF(k) = ( aHybSigmC(k) + aHybSigmC(k-1) )*0.5 _d 0
299     ENDDO
300     ENDIF
301     IF ( .NOT.setCenterDr ) THEN
302     C- Center position at the middle between 2 interfaces
303     DO k=1,Nr
304     bHybSigmC(k) = ( bHybSigmF(k) + bHybSigmF(k+1) )*0.5 _d 0
305     aHybSigmC(k) = ( aHybSigmF(k) + aHybSigmF(k+1) )*0.5 _d 0
306     ENDDO
307     ENDIF
308    
309     C-- Vertical increment:
310     DO k=1,Nr
311     dAHybSigF(k) = ( aHybSigmF(k+1) - aHybSigmF(k) )*rkSign
312     dBHybSigF(k) = ( bHybSigmF(k+1) - bHybSigmF(k) )*rkSign
313     ENDDO
314     DO k=2,Nr
315     dAHybSigC(k) = ( aHybSigmC(k) - aHybSigmC(k-1) )*rkSign
316     dBHybSigC(k) = ( bHybSigmC(k) - bHybSigmC(k-1) )*rkSign
317     ENDDO
318     dAHybSigC(1) = ( aHybSigmC(1) - aHybSigmF(1) )*rkSign
319     dBHybSigC(1) = ( bHybSigmC(1) - bHybSigmF(1) )*rkSign
320     dAHybSigC(Nr+1) = ( aHybSigmF(Nr+1) - aHybSigmC(Nr) )*rkSign
321     dBHybSigC(Nr+1) = ( bHybSigmF(Nr+1) - bHybSigmC(Nr) )*rkSign
322    
323     C-- Check for miss-match between vertical discretization :
324     maxErrC = 0.
325     maxErrF = 0.
326     epsil = 1. _d -9
327     DO k=1,Nr
328     tmpError = ( rC(k)-rF(Nr+1) )*recip_fullDepth
329     & - ( aHybSigmC(k)+bHybSigmC(k) )
330     IF ( ABS(tmpError).GT.epsil ) THEN
331     IF ( maxErrC.LE.epsil ) THEN
332     WRITE(msgBuf,'(2A)') 'S/R INI_VERTICAL_GRID:',
333     & ' rC and Hybrid-Sigma Coeff miss-match'
334     CALL PRINT_ERROR( msgBuf, myThid )
335     ENDIF
336     WRITE(msgBuf,'(A,I4,2(A,1PE14.6),A,1P2E14.6)')
337     & ' k=', k,' , err=', tmpError, ' ; rC=', rC(k),
338     & ' , a & b=', aHybSigmC(k), bHybSigmC(k)
339     CALL PRINT_ERROR( msgBuf, myThid )
340     ENDIF
341     maxErrC = MAX( maxErrC, ABS(tmpError) )
342     ENDDO
343     DO k=1,Nr+1
344     tmpError = ( rF(k)-rF(Nr+1) )*recip_fullDepth
345     & - ( aHybSigmF(k)+bHybSigmF(k) )
346     IF ( ABS(tmpError).GT.epsil ) THEN
347     IF ( maxErrF.LE.epsil ) THEN
348     WRITE(msgBuf,'(2A)') 'S/R INI_VERTICAL_GRID:',
349     & ' rF and Hybrid-Sigma Coeff miss-match'
350     CALL PRINT_ERROR( msgBuf, myThid )
351     ENDIF
352     WRITE(msgBuf,'(A,I4,2(A,1PE14.6),A,1P2E14.6)')
353     & ' k=', k,' , err=', tmpError, ' ; rF=', rF(k),
354     & ' , a & b=', aHybSigmF(k), bHybSigmF(k)
355     CALL PRINT_ERROR( msgBuf, myThid )
356     ENDIF
357     maxErrF = MAX( maxErrF, ABS(tmpError) )
358     ENDDO
359     WRITE(msgBuf,'(2A,1PE14.6)') 'S/R INI_VERTICAL_GRID:',
360     & ' matching of aHybSigmC & rC :', maxErrC
361     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
362     & SQUEEZE_RIGHT, myThid )
363     WRITE(msgBuf,'(2A,1PE14.6)') 'S/R INI_VERTICAL_GRID:',
364     & ' matching of aHybSigmF & rF :', maxErrF
365     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
366     & SQUEEZE_RIGHT, myThid )
367     IF ( maxErrC.GT.epsil .OR. maxErrF.GT.epsil ) THEN
368     WRITE(msgBuf,'(2A)') 'S/R INI_VERTICAL_GRID:',
369     & ' rC,rF and Hybrid-Sigma Coeff miss-match'
370     CALL PRINT_ERROR( msgBuf, myThid )
371     STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID'
372     ENDIF
373    
374     C--- End setting-up Hybrid-Sigma vertical coordinate:
375     ENDIF
376    
377 jmc 1.15 _END_MASTER(myThid)
378     _BARRIER
379    
380 adcroft 1.1 RETURN
381     END

  ViewVC Help
Powered by ViewVC 1.1.22