/[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.17 by jmc, Wed Nov 29 04:39:06 2006 UTC revision 1.21 by jmc, Thu Aug 18 20:29:49 2011 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)
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.  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  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.  C     coordinate (pressure or meters) relative to the vertical index (k) is important.
# Line 50  C     rkSign =  1 applies when k and the Line 59  C     rkSign =  1 applies when k and the
59           gravitySign = 1. _d 0           gravitySign = 1. _d 0
60        ENDIF        ENDIF
61    
62        IF ( .NOT.(setCenterDr.OR.setInterFDr) ) THEN        IF ( .NOT.(setInterFDr.OR.setCenterDr) ) THEN
63           WRITE(msgBuf,'(A)')          WRITE(msgBuf,'(A)')
64       &  'S/R INI_VERTICAL_GRID: neither delR nor delRc are defined'       &  'S/R INI_VERTICAL_GRID: neither delR nor delRc are defined'
65           CALL PRINT_ERROR( msgBuf, myThid )          CALL PRINT_ERROR( msgBuf, myThid )
66           WRITE(msgBuf,'(A)')          WRITE(msgBuf,'(A)')
67       &  'S/R INI_VERTICAL_GRID: Need at least 1 of the 2 (delR,delRc)'       &  'S/R INI_VERTICAL_GRID: Need at least 1 of the 2 (delR,delRc)'
68           CALL PRINT_ERROR( msgBuf, myThid )          CALL PRINT_ERROR( msgBuf, myThid )
69           STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID'          STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID'
70        ENDIF        ENDIF
71    
72  C---  Set Level r-thickness (drF) and Center r-distances (drC)  C---  Set Level r-thickness (drF) and Center r-distances (drC)
73    
74        IF (setInterFDr) THEN        IF (setInterFDr) THEN
75  C--   Interface r-distances are defined:  C--   Interface r-distances are defined:
76         DO k=1,Nr         DO k=1,Nr
77           drF(k) = delR(k)           drF(k) = delR(k)
78         ENDDO         ENDDO
79  C-    Check that all thickness are > 0 :  C-    Check that all thickness are > 0 :
80         DO k=1,Nr         DO k=1,Nr
81          IF (delR(k).LE.0.) THEN          IF (delR(k).LE.0.) THEN
82           WRITE(msgBuf,'(A,I4,A,E16.8)')           WRITE(msgBuf,'(A,I4,A,E16.8)')
83       &  'S/R INI_VERTICAL_GRID: delR(k=',k,' )=',delR(k)       &  'S/R INI_VERTICAL_GRID: delR(k=',k,' )=',delR(k)
84           CALL PRINT_ERROR( msgBuf, myThid )           CALL PRINT_ERROR( msgBuf, myThid )
85           WRITE(msgBuf,'(A)')           WRITE(msgBuf,'(A)')
86       &  'S/R INI_VERTICAL_GRID: Vert. grid spacing MUST BE > 0'       &  'S/R INI_VERTICAL_GRID: Vert. grid spacing MUST BE > 0'
87           CALL PRINT_ERROR( msgBuf, myThid )           CALL PRINT_ERROR( msgBuf, myThid )
88           STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID'           STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID'
89          ENDIF          ENDIF
90         ENDDO         ENDDO
91        ELSE        ELSE
92  C--   Interface r-distances undefined:  C--   Interface r-distances undefined:
93  C     assume Interface at middle between 2 Center  C     assume Interface at middle between 2 Center
94         drF(1) = delRc(1)         drF(1) = delRc(1)
95         DO k=2,Nr         DO k=2,Nr
96           drF(k-1) = 0.5 _d 0 *delRc(k) + drF(k-1)  c        drF(k-1) = 0.5 _d 0 *delRc(k) + drF(k-1)
97           drF( k ) = 0.5 _d 0 *delRc(k)  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         ENDDO
103         drF(Nr) = delRc(Nr+1) + drF(Nr)         drF(Nr) = delRc(Nr+1) + drF(Nr)
104        ENDIF        ENDIF
# Line 155  C-    Calculate reciprol vertical grid s Line 169  C-    Calculate reciprol vertical grid s
169         recip_drF(k)   = 1. _d 0/drF(k)         recip_drF(k)   = 1. _d 0/drF(k)
170        ENDDO        ENDDO
171    
172    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        _END_MASTER(myThid)        _END_MASTER(myThid)
382        _BARRIER        _BARRIER
383    
Legend:
Removed from v.1.17  
changed lines
  Added in v.1.21
  ViewVC Help
Powered by ViewVC 1.1.22