/[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.16 by jmc, Tue Nov 28 22:49:29 2006 UTC revision 1.17 by jmc, Wed Nov 29 04:39:06 2006 UTC
# Line 7  CBOP Line 7  CBOP
7  C     !ROUTINE: INI_VERTICAL_GRID  C     !ROUTINE: INI_VERTICAL_GRID
8  C     !INTERFACE:  C     !INTERFACE:
9        SUBROUTINE INI_VERTICAL_GRID( myThid )        SUBROUTINE INI_VERTICAL_GRID( myThid )
10    
11  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
12  C     *==========================================================*  C     *==========================================================*
13  C     | SUBROUTINE INI_VERTICAL_GRID                                C     | SUBROUTINE INI_VERTICAL_GRID
14  C     | o Initialise vertical gridding arrays                      C     | o Initialise vertical gridding arrays
15  C     *==========================================================*  C     *==========================================================*
16  C     \ev  C     \ev
17    
# Line 24  C     === Global variables === Line 25  C     === Global variables ===
25    
26  C     !INPUT/OUTPUT PARAMETERS:  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
30    
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 meesage buffer
35        INTEGER K        INTEGER k
36          _RL     tmpRatio, checkRatio1, checkRatio2
37        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
38  CEOP  CEOP
39    
# Line 39  CEOP Line 41  CEOP
41    
42  C--   Set factors required for mixing pressure and meters as vertical coordinate.  C--   Set factors required for mixing pressure and meters as vertical coordinate.
43  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
44  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.
45  C     rkSign = -1 applies when K and the coordinate are in the opposite sense.  C     rkSign = -1 applies when k and the coordinate are in the opposite sense.
46  C     rkSign =  1 applies when K and the coordinate are in the same sense.  C     rkSign =  1 applies when k and the coordinate are in the same sense.
47        rkSign       = -1. _d 0        rkSign       = -1. _d 0
48        gravitySign  = -1. _d 0        gravitySign  = -1. _d 0
49        IF ( usingPCoords ) THEN        IF ( usingPCoords ) THEN
50           gravitySign = 1. _d 0           gravitySign = 1. _d 0
51        ENDIF        ENDIF
52    
53        IF (setCenterDr) THEN        IF ( .NOT.(setCenterDr.OR.setInterFDr) ) THEN
54  C-- Interface at middle between 2 centers :           WRITE(msgBuf,'(A)')
55         &  'S/R INI_VERTICAL_GRID: neither delR nor delRc are defined'
56             CALL PRINT_ERROR( msgBuf, myThid )
57             WRITE(msgBuf,'(A)')
58         &  'S/R INI_VERTICAL_GRID: Need at least 1 of the 2 (delR,delRc)'
59             CALL PRINT_ERROR( msgBuf, myThid )
60             STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID'
61          ENDIF
62    
63    C---  Set Level r-thickness (drF) and Center r-distances (drC)
64          IF (setInterFDr) THEN
65    C--   Interface r-distances are defined:
66           DO k=1,Nr
67             drF(k) = delR(k)
68           ENDDO
69  C-    Check that all thickness are > 0 :  C-    Check that all thickness are > 0 :
70         DO K=1,Nr+1         DO k=1,Nr
71          IF (delRc(K).LE.0.) THEN          IF (delR(k).LE.0.) THEN
72           WRITE(msgBuf,'(A,I4,A,E16.8)')           WRITE(msgBuf,'(A,I4,A,E16.8)')
73       &  'S/R INI_VERTICAL_GRID: delRc(K=',K,' )=',delRc(K)       &  'S/R INI_VERTICAL_GRID: delR(k=',k,' )=',delR(k)
74           CALL PRINT_ERROR( msgBuf , 1)           CALL PRINT_ERROR( msgBuf, myThid )
75           WRITE(msgBuf,'(A)')           WRITE(msgBuf,'(A)')
76       &  'S/R INI_VERTICAL_GRID: Vert. grid spacing MUST BE > 0'       &  'S/R INI_VERTICAL_GRID: Vert. grid spacing MUST BE > 0'
77           CALL PRINT_ERROR( msgBuf , 1)           CALL PRINT_ERROR( msgBuf, myThid )
78           STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID'           STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID'
79          ENDIF          ENDIF
        ENDDO  
   
 C-    Calculate depths of centers and interfaces  
         rF(1)  = Ro_SeaLevel  
         rC(1)  = rF(1) + rkSign*delRc(1)  
         drC(1) = delRc(1)  
         drF(1) = delRc(1)  
        DO K=2,Nr  
         drC(K)   = delRc(K)  
         drF(K-1) =  drF(K-1) + 0.5 _d 0*delRc(K)  
         drF(K)   = 0.5 _d 0*delRc(K)  
         rC(K)    = rC(K-1) + rkSign*drC(K)  
         rF(K)    = rF(K-1) + rkSign*drF(K-1)  
80         ENDDO         ENDDO
         drF(Nr)  = drF(Nr) + delRc(Nr+1)  
         rF(Nr+1) = rF(Nr) + rkSign*drF(Nr)  
   
81        ELSE        ELSE
82  C-- Center at middle between 2 interfaces :  C--   Interface r-distances undefined:
83    C     assume Interface at middle between 2 Center
84           drF(1) = delRc(1)
85           DO k=2,Nr
86             drF(k-1) = 0.5 _d 0 *delRc(k) + drF(k-1)
87             drF( k ) = 0.5 _d 0 *delRc(k)
88           ENDDO
89           drF(Nr) = delRc(Nr+1) + drF(Nr)
90          ENDIF
91    
92          IF (setCenterDr) THEN
93    C--   Cell Center r-distances are defined:
94           DO k=1,Nr
95             drC(k) = delRc(k)
96           ENDDO
97  C-    Check that all thickness are > 0 :  C-    Check that all thickness are > 0 :
98         DO K=1,Nr         DO k=1,Nr+1
99          IF (delR(K).LE.0.) THEN          IF (delRc(k).LE.0.) THEN
100           WRITE(msgBuf,'(A,I4,A,E16.8)')           WRITE(msgBuf,'(A,I4,A,E16.8)')
101       &  'S/R INI_VERTICAL_GRID: delR(K=',K,' )=',delR(K)       &  'S/R INI_VERTICAL_GRID: delRc(k=',k,' )=',delRc(k)
102           CALL PRINT_ERROR( msgBuf , 1)           CALL PRINT_ERROR( msgBuf, myThid )
103           WRITE(msgBuf,'(A)')           WRITE(msgBuf,'(A)')
104       &  'S/R INI_VERTICAL_GRID: Vert. grid spacing MUST BE > 0'       &  'S/R INI_VERTICAL_GRID: Vert. grid spacing MUST BE > 0'
105           CALL PRINT_ERROR( msgBuf , 1)           CALL PRINT_ERROR( msgBuf, myThid )
106           STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID'           STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID'
107          ENDIF          ENDIF
108         ENDDO         ENDDO
109          ELSE
110    C--   Cell Center r-distances undefined:
111    C     assume Center at middle between 2 Interfaces
112           drC(1)  = 0.5 _d 0 *delR(1)
113           DO k=2,Nr
114             drC(k) = 0.5 _d 0 *(delR(k-1)+delR(k))
115           ENDDO
116          ENDIF
117    
118  C-    Calculate depths of interfaces and centers  C---  Set r-position of  interFace (rF) and cell-Center (rC):
119        rF(1) = Ro_SeaLevel        rF(1)    = Ro_SeaLevel
120        DO K=1,Nr        DO k=1,Nr
121         drF(K)     = delR(K)         rF(k+1) = rF(k)  + rkSign*drF(k)
        rF(K+1) = rF(K) + rkSign*delR(K)  
122        ENDDO        ENDDO
123        drC(1)      = delR(1) * 0.5 _d 0        rC(1)   = rF(1)   + rkSign*drC(1)
124        rC(1)       = rf(1) + rkSign*delR(1) * 0.5 _d 0        DO k=2,Nr
125        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) + rkSign*drC(K)  
126        ENDDO        ENDDO
127    
128  C--  C---  Check vertical discretization :
129        ENDIF        checkRatio2 = 100.
130          checkRatio1 = 1. _d 0 / checkRatio2
131          DO k=1,Nr
132           tmpRatio = 0.
133           IF ( (rC(k)-rF(k+1)) .NE. 0. )
134         &   tmpRatio = (rF(k)-rC(k)) / (rC(k)-rF(k+1))
135           IF ( tmpRatio.LT.checkRatio1 .OR. tmpRatio.GT.checkRatio2 ) THEN
136    c       write(0,*) 'drF=',drF
137    c       write(0,*) 'drC=',drC
138    c       write(0,*) 'rF=',rF
139    c       write(0,*) 'rC=',rC
140            WRITE(msgBuf,'(A,I4,A,E16.8)')
141         &   'S/R INI_VERTICAL_GRID: Invalid relative position, level k=',
142         &     k, ' :', tmpRatio
143            CALL PRINT_ERROR( msgBuf, myThid )
144            WRITE(msgBuf,'(A,1PE14.6,A,2E14.6)')
145         &   'S/R INI_VERTICAL_GRID: rC=', rC(k),
146         &   ' , rF(k,k+1)=',rF(k),rF(k+1)
147            CALL PRINT_ERROR( msgBuf, myThid )
148            STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID'
149           ENDIF
150          ENDDO
151    
152  C-    Calculate reciprol vertical grid spacing :  C-    Calculate reciprol vertical grid spacing :
153        DO K=1,Nr        DO k=1,Nr
154         recip_drC(K)   = 1. _d 0/drC(K)         recip_drC(k)   = 1. _d 0/drC(k)
155         recip_drF(K)   = 1. _d 0/drF(K)         recip_drF(k)   = 1. _d 0/drF(k)
156        ENDDO        ENDDO
157    
158        _END_MASTER(myThid)        _END_MASTER(myThid)
Legend:
Removed from v.1.16  
changed lines
  Added in v.1.17
  ViewVC Help
Powered by ViewVC 1.1.22