model/src/ini_pressure.F

C $Header: 
C $Name:  

#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: INI_PRESSURE
C     !INTERFACE:
      SUBROUTINE INI_PRESSURE( myThid )
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE INI_PRESSURE
C     | o initialise the pressure field consistently with 
C     |   temperature and salinity
C     |   - needs to be called after ini_theta, ini_salt, and 
C     |     ini_psurf
C     |   - does not include surface pressure loading, because
C     |     that is only available until after 
C     |     CALL packages_init_variables
C     *==========================================================*
C     \ev

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

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     myThid -  Number of this instance of INI_PRESSURE
      INTEGER myThid

C     !LOCAL VARIABLES:
C     == Local variables ==
C     dPhiHydX,Y :: Gradient (X & Y directions) of Hyd. Potential
C     bi,bj  - Loop counters
C     I,J,K
      INTEGER bi, bj
      INTEGER  I,  J, K
      INTEGER ncount, iMin, iMax, jMin, jMax, npiter
      _RL phiHyd(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr) 
      _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL pold(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
      _RL rmspp, rmsppold 

      CHARACTER*(MAX_LEN_MBUF) msgBuf
CEOP
      iMin = 1-OLx
      iMax = sNx+OLx
      jMin = 1-OLy
      jMax = sNy+OLy
      
      WRITE(msgBuf,'(a)')
     &     'Start initial hydrostatic pressure computation'
      CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, 
     &     SQUEEZE_RIGHT , 1)
C     initialise pressure to water column thickness times a constant
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO k = 1,Nr
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           phiHyd(i,j,k) = 0. _d 0
           dPhiHydX(i,j) = 0. _d 0
           dPhiHydY(i,j) = 0. _d 0
          ENDDO
         ENDDO
         CALL STORE_PRESSURE(  bi, bj, k, phiHyd, myThid )
        ENDDO
       ENDDO
      ENDDO

      IF ( buoyancyRelation .EQ. 'OCEANIC' ) THEN

       rmspp    = 1. _d 0
       rmsppold = 0. _d 0
       npiter = 0

C     iterate pressure p = integral of (g*rho(p)*dz)
       DO npiter = 1, 15
        rmsppold = rmspp
        rmspp    = 0. _d 0
        ncount   = 0
        DO bj = myByLo(myThid), myByHi(myThid)
         DO bi = myBxLo(myThid), myBxHi(myThid)
          DO k = 1, Nr
C     for each level save old pressure and compute new pressure 
           DO j=jMin,jMax
            DO i=iMin,iMax
             pold(i,j,k) = pressure(i,j,k,bi,bj)
            ENDDO
           ENDDO
           CALL CALC_PHI_HYD(
     I          bi, bj, iMin, iMax, jMin, jMax, k,
     I          theta, salt, 
     U          phiHyd, 
     O          dPhiHydX, dPhiHydY,
     I          startTime, nIter0, myThid)
           CALL STORE_PRESSURE(  bi, bj, k, phiHyd, myThid )
          ENDDO
C     compute convergence criterion
          DO k = 1, Nr
           DO j=jMin,jMax
            DO i=iMin,iMax
             rmspp = rmspp 
     &            + (pressure(i,j,k,bi,bj)-pold(i,j,k))**2
             ncount = ncount + int(maskC(i,j,k,bi,bj))
            ENDDO
           ENDDO
          ENDDO
         ENDDO
        ENDDO
Cml        WRITE(msgBuf,'(I10.10)') npiter
Cml        CALL WRITE_FLD_XYZ_RL( 'POLD.',msgBuf,pold,npiter,myThid)
Cml        CALL WRITE_FLD_XYZ_RL( 'PINI.',msgBuf,pressure,npiter,myThid)
        IF ( ncount .EQ. 0 ) THEN
           rmspp = 0. _d 0
        ELSE
           rmspp = sqrt(rmspp/ncount)
        ENDIF
        WRITE(msgBuf,'(a,i2,a,e20.13)')
     &       'Iteration ', npiter, ', RMS-difference = ', rmspp
        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, 
     &       SQUEEZE_RIGHT , 1)

       ENDDO
C     print some diagnostics    
       IF ( rmspp .ne. 0. ) THEN
        IF ( rmspp .EQ. rmsppold ) THEN
         WRITE(msgBuf,'(a)')
     &      'Initial hydrostatic pressure did not converge perfectly,' 
         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, 
     &        SQUEEZE_RIGHT , 1)
         WRITE(msgBuf,'(a)')
     &      'but the RMS-difference is constant, indicating that the'
         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, 
     &        SQUEEZE_RIGHT , 1)
         WRITE(msgBuf,'(a)')
     &      'algorithm converged within machine precision.'    
         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, 
     &        SQUEEZE_RIGHT , 1)
        ELSE
         WRITE(msgBuf,'(a,i2,a)')
     &       'Initial hydrostatic pressure did not converge after ', 
     &          npiter-1, ' steps'
         CALL PRINT_ERROR( msgBuf, myThid )
         STOP 'ABNORMAL END: S/R INI_PRESSURE'
        ENDIF
       ENDIF
       WRITE(msgBuf,'(a)')
     &     'Initial hydrostatic pressure converged.' 
       CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, 
     &      SQUEEZE_RIGHT , 1)

      ELSE
C     print a message and DO nothing
       WRITE(msgBuf,'(a,a)')
     &        'Pressure is predetermined for buoyancyRelation ',
     &        buoyancyRelation(1:11)
       CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, 
     &      SQUEEZE_RIGHT , 1)

      ENDIF

      WRITE(msgBuf,'(a)') ' '
      CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, 
     &     SQUEEZE_RIGHT , 1)

      RETURN 
      END
1	mlosch	1.1	C $Header:
2			C $Name:
3
4			#include "CPP_OPTIONS.h"
5
6			CBOP
7			C !ROUTINE: INI_PRESSURE
8			C !INTERFACE:
9	jmc	1.3	SUBROUTINE INI_PRESSURE( myThid )
10	mlosch	1.1	C !DESCRIPTION: \bv
11			C ==========================================================
12			C \| SUBROUTINE INI_PRESSURE
13			C \| o initialise the pressure field consistently with
14			C \| temperature and salinity
15			C \| - needs to be called after ini_theta, ini_salt, and
16			C \| ini_psurf
17			C \| - does not include surface pressure loading, because
18			C \| that is only available until after
19	jmc	1.3	C \| CALL packages_init_variables
20	mlosch	1.1	C ==========================================================
21			C \ev
22
23			C !USES:
24	jmc	1.3	IMPLICIT NONE
25	mlosch	1.1	C == Global variables ==
26			#include "SIZE.h"
27			#include "EEPARAMS.h"
28			#include "PARAMS.h"
29			#include "EOS.h"
30			#include "GRID.h"
31			#include "DYNVARS.h"
32
33			C !INPUT/OUTPUT PARAMETERS:
34			C == Routine arguments ==
35			C myThid - Number of this instance of INI_PRESSURE
36			INTEGER myThid
37
38			C !LOCAL VARIABLES:
39			C == Local variables ==
40	jmc	1.3	C dPhiHydX,Y :: Gradient (X & Y directions) of Hyd. Potential
41	mlosch	1.1	C bi,bj - Loop counters
42			C I,J,K
43			INTEGER bi, bj
44			INTEGER I, J, K
45			INTEGER ncount, iMin, iMax, jMin, jMax, npiter
46	jmc	1.3	_RL phiHyd(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
47			_RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
48			_RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
49	mlosch	1.1	_RL pold(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
50			_RL rmspp, rmsppold
51
52			CHARACTER*(MAX_LEN_MBUF) msgBuf
53	mlosch	1.2	CEOP
54	mlosch	1.1	iMin = 1-OLx
55			iMax = sNx+OLx
56			jMin = 1-OLy
57			jMax = sNy+OLy
58
59	jmc	1.3	WRITE(msgBuf,'(a)')
60	mlosch	1.1	& 'Start initial hydrostatic pressure computation'
61			CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
62			& SQUEEZE_RIGHT , 1)
63			C initialise pressure to water column thickness times a constant
64	jmc	1.3	DO bj = myByLo(myThid), myByHi(myThid)
65			DO bi = myBxLo(myThid), myBxHi(myThid)
66			DO k = 1,Nr
67			DO j=1-OLy,sNy+OLy
68			DO i=1-OLx,sNx+OLx
69	mlosch	1.1	phiHyd(i,j,k) = 0. _d 0
70	jmc	1.3	dPhiHydX(i,j) = 0. _d 0
71			dPhiHydY(i,j) = 0. _d 0
72			ENDDO
73			ENDDO
74			CALL STORE_PRESSURE( bi, bj, k, phiHyd, myThid )
75			ENDDO
76			ENDDO
77			ENDDO
78	mlosch	1.1
79	jmc	1.3	IF ( buoyancyRelation .EQ. 'OCEANIC' ) THEN
80	mlosch	1.1
81			rmspp = 1. _d 0
82			rmsppold = 0. _d 0
83			npiter = 0
84
85			C iterate pressure p = integral of (grho(p)dz)
86	jmc	1.3	DO npiter = 1, 15
87	mlosch	1.1	rmsppold = rmspp
88			rmspp = 0. _d 0
89			ncount = 0
90	jmc	1.3	DO bj = myByLo(myThid), myByHi(myThid)
91			DO bi = myBxLo(myThid), myBxHi(myThid)
92			DO k = 1, Nr
93	mlosch	1.1	C for each level save old pressure and compute new pressure
94	jmc	1.3	DO j=jMin,jMax
95			DO i=iMin,iMax
96	mlosch	1.1	pold(i,j,k) = pressure(i,j,k,bi,bj)
97	jmc	1.3	ENDDO
98			ENDDO
99			CALL CALC_PHI_HYD(
100			I bi, bj, iMin, iMax, jMin, jMax, k,
101			I theta, salt,
102			U phiHyd,
103			O dPhiHydX, dPhiHydY,
104			I startTime, nIter0, myThid)
105			CALL STORE_PRESSURE( bi, bj, k, phiHyd, myThid )
106			ENDDO
107	mlosch	1.1	C compute convergence criterion
108	jmc	1.3	DO k = 1, Nr
109			DO j=jMin,jMax
110			DO i=iMin,iMax
111	mlosch	1.1	rmspp = rmspp
112			& + (pressure(i,j,k,bi,bj)-pold(i,j,k))**2
113			ncount = ncount + int(maskC(i,j,k,bi,bj))
114	jmc	1.3	ENDDO
115			ENDDO
116			ENDDO
117			ENDDO
118			ENDDO
119	mlosch	1.1	Cml WRITE(msgBuf,'(I10.10)') npiter
120			Cml CALL WRITE_FLD_XYZ_RL( 'POLD.',msgBuf,pold,npiter,myThid)
121			Cml CALL WRITE_FLD_XYZ_RL( 'PINI.',msgBuf,pressure,npiter,myThid)
122	jmc	1.3	IF ( ncount .EQ. 0 ) THEN
123	mlosch	1.1	rmspp = 0. _d 0
124	jmc	1.3	ELSE
125	mlosch	1.1	rmspp = sqrt(rmspp/ncount)
126	jmc	1.3	ENDIF
127			WRITE(msgBuf,'(a,i2,a,e20.13)')
128	mlosch	1.1	& 'Iteration ', npiter, ', RMS-difference = ', rmspp
129			CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
130			& SQUEEZE_RIGHT , 1)
131
132	jmc	1.3	ENDDO
133	mlosch	1.1	C print some diagnostics
134	jmc	1.3	IF ( rmspp .ne. 0. ) THEN
135			IF ( rmspp .EQ. rmsppold ) THEN
136			WRITE(msgBuf,'(a)')
137	mlosch	1.1	& 'Initial hydrostatic pressure did not converge perfectly,'
138			CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
139			& SQUEEZE_RIGHT , 1)
140	jmc	1.3	WRITE(msgBuf,'(a)')
141	mlosch	1.1	& 'but the RMS-difference is constant, indicating that the'
142			CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
143			& SQUEEZE_RIGHT , 1)
144	jmc	1.3	WRITE(msgBuf,'(a)')
145	mlosch	1.1	& 'algorithm converged within machine precision.'
146			CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
147			& SQUEEZE_RIGHT , 1)
148	jmc	1.3	ELSE
149			WRITE(msgBuf,'(a,i2,a)')
150	mlosch	1.1	& 'Initial hydrostatic pressure did not converge after ',
151			& npiter-1, ' steps'
152			CALL PRINT_ERROR( msgBuf, myThid )
153			STOP 'ABNORMAL END: S/R INI_PRESSURE'
154	jmc	1.3	ENDIF
155			ENDIF
156			WRITE(msgBuf,'(a)')
157	mlosch	1.1	& 'Initial hydrostatic pressure converged.'
158			CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
159			& SQUEEZE_RIGHT , 1)
160
161	jmc	1.3	ELSE
162			C print a message and DO nothing
163			WRITE(msgBuf,'(a,a)')
164	mlosch	1.1	& 'Pressure is predetermined for buoyancyRelation ',
165			& buoyancyRelation(1:11)
166			CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
167			& SQUEEZE_RIGHT , 1)
168
169	jmc	1.3	ENDIF
170	mlosch	1.1
171	jmc	1.3	WRITE(msgBuf,'(a)') ' '
172	mlosch	1.1	CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
173			& SQUEEZE_RIGHT , 1)
174
175	jmc	1.3	RETURN
176			END