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

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

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

revision 1.9 by adcroft, Tue May 18 17:42:35 1999 UTC revision 1.22 by jmc, Fri Jul 26 14:22:22 2013 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    
6  CStartOfInterface  CBOP
7    C     !ROUTINE: INI_VERTICAL_GRID
8    C     !INTERFACE:
9        SUBROUTINE INI_VERTICAL_GRID( myThid )        SUBROUTINE INI_VERTICAL_GRID( myThid )
 C     /==========================================================\  
 C     | SUBROUTINE INI_VERTICAL_GRID                             |  
 C     | o Initialise vertical gridding arrays                    |  
 C     |==========================================================|  
 C     |                                                          |  
 C     \==========================================================/  
       IMPLICIT NONE  
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 ===  C     === Global variables ===
21  #include "SIZE.h"  #include "SIZE.h"
22  #include "EEPARAMS.h"  #include "EEPARAMS.h"
23  #include "PARAMS.h"  #include "PARAMS.h"
24  #include "GRID.h"  #include "GRID.h"
25    
26    C     !INPUT/OUTPUT PARAMETERS:
27  C     == Routine arguments ==  C     == Routine arguments ==
28  C     myThid -  Number of this instance of INI_DEPTHS  C     myThid   :: my Thread Id number
29        INTEGER myThid        INTEGER myThid
 CEndOfInterface  
30    
31    C     !LOCAL VARIABLES:
32  C     == Local variables ==  C     == Local variables ==
33  C     K  C     k        :: loop index
34        INTEGER K  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  C     Calculate depths of centers and interfaces        IF ( .NOT.(setInterFDr.OR.setCenterDr) ) THEN
63        rF(1) = 0. _d 0          WRITE(msgBuf,'(A)')
64        DO K=1,Nr       &  'S/R INI_VERTICAL_GRID: neither delR nor delRc are defined'
65         drF(K)     = delR(K)          CALL PRINT_ERROR( msgBuf, myThid )
66         rF(K+1) = rF(K)-rkFac*delR(K)          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        ENDDO
138        drC(1)      = delR(1) * 0.5 _d 0        rC(1)   = rF(1)   + rkSign*drC(1)
139        rC(1)       = rf(1)-rkFac*delR(1) * 0.5 _d 0        DO k=2,Nr
140        DO K=2,Nr          rC(k) = rC(k-1) + rkSign*drC(k)
        drC(K)     = 0.5 _d 0 *(delR(K-1)+delR(K))  
        rC(K)      = rC(K-1) - rkFac*drC(K)  
141        ENDDO        ENDDO
142        DO K=1,Nr  
143         saFac(K)  = 1. _d 0  C---  Check vertical discretization :
144         recip_drC(K)   = 1. _d 0/drC(K)        checkRatio2 = 100.
145         recip_drF(K)   = 1. _d 0/drF(K)        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        ENDDO
166  C  
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        RETURN
388        END        END
Legend:
Removed from v.1.9  
changed lines
  Added in v.1.22
  ViewVC Help
Powered by ViewVC 1.1.22