model/src/ini_vertical_grid.F

C $Header: /u/gcmpack/MITgcm/model/src/ini_vertical_grid.F,v 1.18 2006/12/05 05:25:08 jmc Exp $
C $Name:  $

#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: INI_VERTICAL_GRID
C     !INTERFACE:
      SUBROUTINE INI_VERTICAL_GRID( myThid )

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE INI_VERTICAL_GRID
C     | o Initialise vertical gridding arrays
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     myThid   :: my Thread Id number
      INTEGER myThid

C     !LOCAL VARIABLES:
C     == Local variables ==
C     k        :: loop index
C     msgBuf   :: Informational/error meesage buffer
      INTEGER k
      _RL     tmpRatio, checkRatio1, checkRatio2
      CHARACTER*(MAX_LEN_MBUF) msgBuf
CEOP

      _BEGIN_MASTER(myThid)

      WRITE(msgBuf,'(A,2(A,L5))') 'Enter INI_VERTICAL_GRID:',
     &                            ' setInterFDr=', setInterFDr,
     &                          ' ; setCenterDr=', setCenterDr
      CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                    SQUEEZE_RIGHT, myThid )

C--   Set factors required for mixing pressure and meters as vertical coordinate.
C     rkSign is a "sign" parameter which is used where the orientation of the vertical
C     coordinate (pressure or meters) relative to the vertical index (k) is important.
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 same sense.
      rkSign       = -1. _d 0
      gravitySign  = -1. _d 0
      IF ( usingPCoords ) THEN
         gravitySign = 1. _d 0
      ENDIF

      IF ( .NOT.(setInterFDr.OR.setCenterDr) ) THEN
        WRITE(msgBuf,'(A)')
     &  'S/R INI_VERTICAL_GRID: neither delR nor delRc are defined'
        CALL PRINT_ERROR( msgBuf, myThid )
        WRITE(msgBuf,'(A)')
     &  'S/R INI_VERTICAL_GRID: Need at least 1 of the 2 (delR,delRc)'
        CALL PRINT_ERROR( msgBuf, myThid )
        STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID'
      ENDIF

C---  Set Level r-thickness (drF) and Center r-distances (drC)

      IF (setInterFDr) THEN
C--   Interface r-distances are defined:
       DO k=1,Nr
         drF(k) = delR(k)
       ENDDO
C-    Check that all thickness are > 0 :
       DO k=1,Nr
        IF (delR(k).LE.0.) THEN
         WRITE(msgBuf,'(A,I4,A,E16.8)')
     &  'S/R INI_VERTICAL_GRID: delR(k=',k,' )=',delR(k)
         CALL PRINT_ERROR( msgBuf, myThid )
         WRITE(msgBuf,'(A)')
     &  'S/R INI_VERTICAL_GRID: Vert. grid spacing MUST BE > 0'
         CALL PRINT_ERROR( msgBuf, myThid )
         STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID'
        ENDIF
       ENDDO
      ELSE
C--   Interface r-distances undefined:
C     assume Interface at middle between 2 Center
       drF(1) = delRc(1)
       DO k=2,Nr
         drF(k-1) = 0.5 _d 0 *delRc(k) + drF(k-1)
         drF( k ) = 0.5 _d 0 *delRc(k)
       ENDDO
       drF(Nr) = delRc(Nr+1) + drF(Nr)
      ENDIF

      IF (setCenterDr) THEN
C--   Cell Center r-distances are defined:
       DO k=1,Nr
         drC(k) = delRc(k)
       ENDDO
C-    Check that all thickness are > 0 :
       DO k=1,Nr+1
        IF (delRc(k).LE.0.) THEN
         WRITE(msgBuf,'(A,I4,A,E16.8)')
     &  'S/R INI_VERTICAL_GRID: delRc(k=',k,' )=',delRc(k)
         CALL PRINT_ERROR( msgBuf, myThid )
         WRITE(msgBuf,'(A)')
     &  'S/R INI_VERTICAL_GRID: Vert. grid spacing MUST BE > 0'
         CALL PRINT_ERROR( msgBuf, myThid )
         STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID'
        ENDIF
       ENDDO
      ELSE
C--   Cell Center r-distances undefined:
C     assume Center at middle between 2 Interfaces
       drC(1)  = 0.5 _d 0 *delR(1)
       DO k=2,Nr
         drC(k) = 0.5 _d 0 *(delR(k-1)+delR(k))
       ENDDO
      ENDIF

C---  Set r-position of  interFace (rF) and cell-Center (rC):
      rF(1)    = Ro_SeaLevel
      DO k=1,Nr
       rF(k+1) = rF(k)  + rkSign*drF(k)
      ENDDO
      rC(1)   = rF(1)   + rkSign*drC(1)
      DO k=2,Nr
        rC(k) = rC(k-1) + rkSign*drC(k)
      ENDDO

C---  Check vertical discretization :
      checkRatio2 = 100.
      checkRatio1 = 1. _d 0 / checkRatio2
      DO k=1,Nr
       tmpRatio = 0.
       IF ( (rC(k)-rF(k+1)) .NE. 0. )
     &   tmpRatio = (rF(k)-rC(k)) / (rC(k)-rF(k+1))
       IF ( tmpRatio.LT.checkRatio1 .OR. tmpRatio.GT.checkRatio2 ) THEN
c       write(0,*) 'drF=',drF
c       write(0,*) 'drC=',drC
c       write(0,*) 'rF=',rF
c       write(0,*) 'rC=',rC
        WRITE(msgBuf,'(A,I4,A,E16.8)')
     &   'S/R INI_VERTICAL_GRID: Invalid relative position, level k=',
     &     k, ' :', tmpRatio
        CALL PRINT_ERROR( msgBuf, myThid )
        WRITE(msgBuf,'(A,1PE14.6,A,2E14.6)')
     &   'S/R INI_VERTICAL_GRID: rC=', rC(k),
     &   ' , rF(k,k+1)=',rF(k),rF(k+1)
        CALL PRINT_ERROR( msgBuf, myThid )
        STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID'
       ENDIF
      ENDDO

C-    Calculate reciprol vertical grid spacing :
      DO k=1,Nr
       recip_drC(k)   = 1. _d 0/drC(k)
       recip_drF(k)   = 1. _d 0/drF(k)
      ENDDO

      _END_MASTER(myThid)
      _BARRIER

      RETURN
      END
1	jmc	1.19	C $Header: /u/gcmpack/MITgcm/model/src/ini_vertical_grid.F,v 1.18 2006/12/05 05:25:08 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			C msgBuf :: Informational/error meesage buffer
35			INTEGER k
36			_RL tmpRatio, checkRatio1, checkRatio2
37	jmc	1.13	CHARACTER*(MAX_LEN_MBUF) msgBuf
38	cnh	1.12	CEOP
39	adcroft	1.1
40	jmc	1.15	_BEGIN_MASTER(myThid)
41
42	jmc	1.18	WRITE(msgBuf,'(A,2(A,L5))') 'Enter INI_VERTICAL_GRID:',
43			& ' setInterFDr=', setInterFDr,
44			& ' ; setCenterDr=', setCenterDr
45			CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
46			& SQUEEZE_RIGHT, myThid )
47
48	jmc	1.16	C-- Set factors required for mixing pressure and meters as vertical coordinate.
49			C rkSign is a "sign" parameter which is used where the orientation of the vertical
50	jmc	1.17	C coordinate (pressure or meters) relative to the vertical index (k) is important.
51			C rkSign = -1 applies when k and the coordinate are in the opposite sense.
52			C rkSign = 1 applies when k and the coordinate are in the same sense.
53	jmc	1.16	rkSign = -1. _d 0
54			gravitySign = -1. _d 0
55			IF ( usingPCoords ) THEN
56			gravitySign = 1. _d 0
57			ENDIF
58
59	jmc	1.18	IF ( .NOT.(setInterFDr.OR.setCenterDr) ) THEN
60			WRITE(msgBuf,'(A)')
61	jmc	1.17	& 'S/R INI_VERTICAL_GRID: neither delR nor delRc are defined'
62	jmc	1.18	CALL PRINT_ERROR( msgBuf, myThid )
63			WRITE(msgBuf,'(A)')
64	jmc	1.17	& 'S/R INI_VERTICAL_GRID: Need at least 1 of the 2 (delR,delRc)'
65	jmc	1.18	CALL PRINT_ERROR( msgBuf, myThid )
66			STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID'
67	jmc	1.17	ENDIF
68	jmc	1.13
69	jmc	1.17	C--- Set Level r-thickness (drF) and Center r-distances (drC)
70	jmc	1.18
71	jmc	1.17	IF (setInterFDr) THEN
72			C-- Interface r-distances are defined:
73			DO k=1,Nr
74	jmc	1.18	drF(k) = delR(k)
75	jmc	1.17	ENDDO
76	jmc	1.13	C- Check that all thickness are > 0 :
77	jmc	1.17	DO k=1,Nr
78	jmc	1.18	IF (delR(k).LE.0.) THEN
79			WRITE(msgBuf,'(A,I4,A,E16.8)')
80	jmc	1.17	& 'S/R INI_VERTICAL_GRID: delR(k=',k,' )=',delR(k)
81	jmc	1.18	CALL PRINT_ERROR( msgBuf, myThid )
82			WRITE(msgBuf,'(A)')
83	jmc	1.13	& 'S/R INI_VERTICAL_GRID: Vert. grid spacing MUST BE > 0'
84	jmc	1.18	CALL PRINT_ERROR( msgBuf, myThid )
85			STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID'
86			ENDIF
87	jmc	1.17	ENDDO
88			ELSE
89			C-- Interface r-distances undefined:
90			C assume Interface at middle between 2 Center
91			drF(1) = delRc(1)
92			DO k=2,Nr
93	jmc	1.18	drF(k-1) = 0.5 _d 0 *delRc(k) + drF(k-1)
94			drF( k ) = 0.5 _d 0 *delRc(k)
95	jmc	1.13	ENDDO
96	jmc	1.17	drF(Nr) = delRc(Nr+1) + drF(Nr)
97			ENDIF
98	jmc	1.13
99	jmc	1.17	IF (setCenterDr) THEN
100			C-- Cell Center r-distances are defined:
101			DO k=1,Nr
102			drC(k) = delRc(k)
103	jmc	1.13	ENDDO
104			C- Check that all thickness are > 0 :
105	jmc	1.17	DO k=1,Nr+1
106			IF (delRc(k).LE.0.) THEN
107	jmc	1.13	WRITE(msgBuf,'(A,I4,A,E16.8)')
108	jmc	1.17	& 'S/R INI_VERTICAL_GRID: delRc(k=',k,' )=',delRc(k)
109			CALL PRINT_ERROR( msgBuf, myThid )
110	jmc	1.13	WRITE(msgBuf,'(A)')
111			& 'S/R INI_VERTICAL_GRID: Vert. grid spacing MUST BE > 0'
112	jmc	1.17	CALL PRINT_ERROR( msgBuf, myThid )
113	jmc	1.13	STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID'
114			ENDIF
115			ENDDO
116	jmc	1.17	ELSE
117			C-- Cell Center r-distances undefined:
118			C assume Center at middle between 2 Interfaces
119			drC(1) = 0.5 _d 0 *delR(1)
120			DO k=2,Nr
121			drC(k) = 0.5 _d 0 *(delR(k-1)+delR(k))
122			ENDDO
123			ENDIF
124	jmc	1.13
125	jmc	1.17	C--- Set r-position of interFace (rF) and cell-Center (rC):
126			rF(1) = Ro_SeaLevel
127			DO k=1,Nr
128			rF(k+1) = rF(k) + rkSign*drF(k)
129	adcroft	1.1	ENDDO
130	jmc	1.17	rC(1) = rF(1) + rkSign*drC(1)
131			DO k=2,Nr
132			rC(k) = rC(k-1) + rkSign*drC(k)
133	adcroft	1.1	ENDDO
134	jmc	1.13
135	jmc	1.17	C--- Check vertical discretization :
136			checkRatio2 = 100.
137			checkRatio1 = 1. _d 0 / checkRatio2
138			DO k=1,Nr
139			tmpRatio = 0.
140			IF ( (rC(k)-rF(k+1)) .NE. 0. )
141			& tmpRatio = (rF(k)-rC(k)) / (rC(k)-rF(k+1))
142			IF ( tmpRatio.LT.checkRatio1 .OR. tmpRatio.GT.checkRatio2 ) THEN
143			c write(0,*) 'drF=',drF
144			c write(0,*) 'drC=',drC
145			c write(0,*) 'rF=',rF
146			c write(0,*) 'rC=',rC
147			WRITE(msgBuf,'(A,I4,A,E16.8)')
148			& 'S/R INI_VERTICAL_GRID: Invalid relative position, level k=',
149			& k, ' :', tmpRatio
150			CALL PRINT_ERROR( msgBuf, myThid )
151			WRITE(msgBuf,'(A,1PE14.6,A,2E14.6)')
152			& 'S/R INI_VERTICAL_GRID: rC=', rC(k),
153			& ' , rF(k,k+1)=',rF(k),rF(k+1)
154			CALL PRINT_ERROR( msgBuf, myThid )
155			STOP 'ABNORMAL END: S/R INI_VERTICAL_GRID'
156			ENDIF
157			ENDDO
158	jmc	1.13
159			C- Calculate reciprol vertical grid spacing :
160	jmc	1.17	DO k=1,Nr
161			recip_drC(k) = 1. _d 0/drC(k)
162			recip_drF(k) = 1. _d 0/drF(k)
163	adcroft	1.1	ENDDO
164	jmc	1.13
165	jmc	1.15	_END_MASTER(myThid)
166			_BARRIER
167
168	adcroft	1.1	RETURN
169			END