/[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.19 by jmc, Fri Sep 5 20:15:28 2008 UTC revision 1.20 by jmc, Sat Sep 11 21:24:52 2010 UTC
# Line 31  C     myThid   :: my Thread Id number Line 31  C     myThid   :: my Thread Id number
31  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
32  C     == Local variables ==  C     == Local variables ==
33  C     k        :: loop index  C     k        :: loop index
34  C     msgBuf   :: Informational/error meesage buffer  C     msgBuf   :: Informational/error message buffer
35        INTEGER k        INTEGER k
36        _RL     tmpRatio, checkRatio1, checkRatio2        _RL     tmpRatio, checkRatio1, checkRatio2
37        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
38          _RL maxErrC, maxErrF, epsil, tmpError
39          _RL rFullDepth, recip_fullDepth
40          _RS rSigBndRS, tmpRS
41  CEOP  CEOP
42    
43        _BEGIN_MASTER(myThid)        _BEGIN_MASTER(myThid)
# Line 162  C-    Calculate reciprol vertical grid s Line 165  C-    Calculate reciprol vertical grid s
165         recip_drF(k)   = 1. _d 0/drF(k)         recip_drF(k)   = 1. _d 0/drF(k)
166        ENDDO        ENDDO
167    
168    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        _END_MASTER(myThid)        _END_MASTER(myThid)
378        _BARRIER        _BARRIER
379    
Legend:
Removed from v.1.19  
changed lines
  Added in v.1.20
  ViewVC Help
Powered by ViewVC 1.1.22