/[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.21 - (hide annotations) (download)
Thu Aug 18 20:29:49 2011 UTC (12 years, 9 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint63g, checkpoint64, checkpoint63p, checkpoint63q, checkpoint63r, checkpoint63s, checkpoint63l, checkpoint63m, checkpoint63n, checkpoint63o, checkpoint63h, checkpoint63i, checkpoint63j, checkpoint63k, checkpoint63d, checkpoint63e, checkpoint63f, checkpoint63b, checkpoint63c, checkpoint64i, checkpoint64h, checkpoint64k, checkpoint64j, checkpoint64a, checkpoint64c, checkpoint64b, checkpoint64e, checkpoint64d, checkpoint64g, checkpoint64f
Changes since 1.20: +6 -2 lines
minor change (identical instructions) to prevent some compilers to produce
 wrong code when trying to optimise this routine.

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

  ViewVC Help
Powered by ViewVC 1.1.22