model/src/ini_pressure.F

C $Header: /u/gcmpack/MITgcm/model/src/ini_pressure.F,v 1.11 2008/09/22 17:55:16 jmc Exp $
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      :: tile indices
C     i,j,k      :: Loop counters
      INTEGER bi, bj
      INTEGER  i,  j, k
      INTEGER  iMin, iMax, jMin, jMax, npiter
      _RL PhiHydF (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL PhiHydC (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL oldPhi  (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL count, rmspp, rmsppold
      _RL sumTile, rmsTile

      CHARACTER*(MAX_LEN_MBUF) msgBuf
CEOP
      iMin = 1-OLx
      iMax = sNx+OLx
      jMin = 1-OLy
      jMax = sNy+OLy

      _BEGIN_MASTER( myThid )
      WRITE(msgBuf,'(a)')
     &     'Start initial hydrostatic pressure computation'
      CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &     SQUEEZE_RIGHT , myThid)
      _END_MASTER( myThid )

      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
           totPhiHyd(i,j,k,bi,bj) = 0. _d 0
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO

      IF ( useDynP_inEos_Zc ) THEN

#ifndef ALLOW_AUTODIFF_TAMC
cph-- deal with this iterative loop for AD once it will
cph-- really be needed;
cph-- would need storing of totPhiHyd for each npiter

       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
        count    = 0. _d 0
        DO bj = myByLo(myThid), myByHi(myThid)
         DO bi = myBxLo(myThid), myBxHi(myThid)
          rmsTile = 0. _d 0
          sumTile = 0. _d 0
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
             phiHydF(i,j) = 0. _d 0
           ENDDO
          ENDDO
          DO k = 1, Nr
C     for each level save old pressure and compute new pressure
           DO j=jMin,jMax
            DO i=iMin,iMax
             oldPhi(i,j) = totPhiHyd(i,j,k,bi,bj)
            ENDDO
           ENDDO
           CALL CALC_PHI_HYD(
     I          bi, bj, iMin, iMax, jMin, jMax, k,
     I          theta, salt,
     U          phiHydF,
     O          phiHydC, dPhiHydX, dPhiHydY,
     I          startTime, -1, myThid )
C     compute convergence criterion
           DO j=1,sNy
            DO i=1,sNx
             rmsTile = rmsTile
     &            + (totPhiHyd(i,j,k,bi,bj)-oldPhi(i,j))**2
     &             * maskC(i,j,k,bi,bj)
             sumTile = sumTile + maskC(i,j,k,bi,bj)
            ENDDO
           ENDDO
C --      end k loop
          ENDDO
          rmspp = rmspp + rmsTile
          count = count + sumTile
         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)
        _GLOBAL_SUM_RL( rmspp, myThid )
        _GLOBAL_SUM_RL( count, myThid )
        IF ( count .EQ. 0. ) THEN
           rmspp = 0. _d 0
        ELSE
           rmspp = SQRT(rmspp/count)
        ENDIF
        _BEGIN_MASTER( myThid )
        WRITE(msgBuf,'(A,I4,A,1PE20.12)')
     &       'Iteration', npiter, ', RMS-difference =', rmspp
        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &       SQUEEZE_RIGHT , myThid)
        _END_MASTER( myThid )

C --   end npiter loop
       ENDDO
C     print some diagnostics    
       _BEGIN_MASTER( myThid )
       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 , myThid)
         WRITE(msgBuf,'(A)')
     &      'but the RMS-difference is constant, indicating that the'
         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &        SQUEEZE_RIGHT , myThid)
         WRITE(msgBuf,'(A)')
     &      'algorithm converged within machine precision.'
         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &        SQUEEZE_RIGHT , myThid)
        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 , myThid)
       _END_MASTER( myThid )

#endif /* ALLOW_AUTODIFF_TAMC */

c-- else of if useDynP_inEos_Zc
      ELSE
C     print a message and DO nothing
       _BEGIN_MASTER( myThid )
       WRITE(msgBuf,'(A,A)')
     &        'Pressure is predetermined for buoyancyRelation ',
     &        buoyancyRelation(1:11)
       CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &      SQUEEZE_RIGHT , myThid)
       _END_MASTER( myThid )

      ENDIF

      _BEGIN_MASTER( myThid )
      WRITE(msgBuf,'(A)') ' '
      CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &     SQUEEZE_RIGHT , myThid)
      _END_MASTER( myThid )

      RETURN
      END
1	jmc	1.12	C $Header: /u/gcmpack/MITgcm/model/src/ini_pressure.F,v 1.11 2008/09/22 17:55:16 jmc Exp $
2	edhill	1.6	C $Name: $
3	mlosch	1.1
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	jmc	1.10	C \| o initialise the pressure field consistently with
14	mlosch	1.1	C \| temperature and salinity
15	jmc	1.10	C \| - needs to be called after ini_theta, ini_salt, and
16	mlosch	1.1	C \| ini_psurf
17			C \| - does not include surface pressure loading, because
18	jmc	1.10	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	jmc	1.10	C bi,bj :: tile indices
42			C i,j,k :: Loop counters
43	mlosch	1.1	INTEGER bi, bj
44	jmc	1.10	INTEGER i, j, k
45	jmc	1.4	INTEGER iMin, iMax, jMin, jMax, npiter
46			_RL PhiHydF (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
47			_RL PhiHydC (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
48	jmc	1.3	_RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
49			_RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
50	jmc	1.4	_RL oldPhi (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
51	jmc	1.8	_RL count, rmspp, rmsppold
52	jmc	1.10	_RL sumTile, rmsTile
53	mlosch	1.1
54			CHARACTER*(MAX_LEN_MBUF) msgBuf
55	mlosch	1.2	CEOP
56	mlosch	1.1	iMin = 1-OLx
57			iMax = sNx+OLx
58			jMin = 1-OLy
59			jMax = sNy+OLy
60	jmc	1.8
61	jmc	1.9	_BEGIN_MASTER( myThid )
62	jmc	1.3	WRITE(msgBuf,'(a)')
63	mlosch	1.1	& 'Start initial hydrostatic pressure computation'
64	jmc	1.8	CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
65			& SQUEEZE_RIGHT , myThid)
66	jmc	1.9	_END_MASTER( myThid )
67
68	jmc	1.3	DO bj = myByLo(myThid), myByHi(myThid)
69			DO bi = myBxLo(myThid), myBxHi(myThid)
70			DO k = 1,Nr
71			DO j=1-OLy,sNy+OLy
72			DO i=1-OLx,sNx+OLx
73	jmc	1.4	totPhiHyd(i,j,k,bi,bj) = 0. _d 0
74	jmc	1.3	ENDDO
75			ENDDO
76			ENDDO
77			ENDDO
78			ENDDO
79	mlosch	1.1
80	jmc	1.4	IF ( useDynP_inEos_Zc ) THEN
81	mlosch	1.1
82	heimbach	1.5	#ifndef ALLOW_AUTODIFF_TAMC
83			cph-- deal with this iterative loop for AD once it will
84			cph-- really be needed;
85			cph-- would need storing of totPhiHyd for each npiter
86
87	mlosch	1.1	rmspp = 1. _d 0
88			rmsppold = 0. _d 0
89			npiter = 0
90
91			C iterate pressure p = integral of (grho(p)dz)
92	jmc	1.4	DO npiter= 1, 15
93	mlosch	1.1	rmsppold = rmspp
94			rmspp = 0. _d 0
95	jmc	1.10	count = 0. _d 0
96	jmc	1.3	DO bj = myByLo(myThid), myByHi(myThid)
97			DO bi = myBxLo(myThid), myBxHi(myThid)
98	jmc	1.10	rmsTile = 0. _d 0
99			sumTile = 0. _d 0
100			DO j=1-OLy,sNy+OLy
101			DO i=1-OLx,sNx+OLx
102	jmc	1.4	phiHydF(i,j) = 0. _d 0
103			ENDDO
104	jmc	1.10	ENDDO
105	jmc	1.3	DO k = 1, Nr
106	jmc	1.10	C for each level save old pressure and compute new pressure
107	jmc	1.3	DO j=jMin,jMax
108			DO i=iMin,iMax
109	jmc	1.4	oldPhi(i,j) = totPhiHyd(i,j,k,bi,bj)
110	jmc	1.3	ENDDO
111			ENDDO
112			CALL CALC_PHI_HYD(
113			I bi, bj, iMin, iMax, jMin, jMax, k,
114	jmc	1.8	I theta, salt,
115	jmc	1.4	U phiHydF,
116			O phiHydC, dPhiHydX, dPhiHydY,
117	jmc	1.11	I startTime, -1, myThid )
118	mlosch	1.1	C compute convergence criterion
119	jmc	1.10	DO j=1,sNy
120			DO i=1,sNx
121			rmsTile = rmsTile
122	jmc	1.4	& + (totPhiHyd(i,j,k,bi,bj)-oldPhi(i,j))**2
123	jmc	1.10	& * maskC(i,j,k,bi,bj)
124			sumTile = sumTile + maskC(i,j,k,bi,bj)
125	jmc	1.3	ENDDO
126			ENDDO
127	jmc	1.4	C -- end k loop
128	jmc	1.3	ENDDO
129	jmc	1.10	rmspp = rmspp + rmsTile
130			count = count + sumTile
131	jmc	1.3	ENDDO
132			ENDDO
133	mlosch	1.1	Cml WRITE(msgBuf,'(I10.10)') npiter
134			Cml CALL WRITE_FLD_XYZ_RL( 'POLD.',msgBuf,pold,npiter,myThid)
135			Cml CALL WRITE_FLD_XYZ_RL( 'PINI.',msgBuf,pressure,npiter,myThid)
136	jmc	1.12	_GLOBAL_SUM_RL( rmspp, myThid )
137			_GLOBAL_SUM_RL( count, myThid )
138	jmc	1.4	IF ( count .EQ. 0. ) THEN
139	mlosch	1.1	rmspp = 0. _d 0
140	jmc	1.3	ELSE
141	jmc	1.10	rmspp = SQRT(rmspp/count)
142	jmc	1.3	ENDIF
143	jmc	1.9	_BEGIN_MASTER( myThid )
144	jmc	1.10	WRITE(msgBuf,'(A,I4,A,1PE20.12)')
145			& 'Iteration', npiter, ', RMS-difference =', rmspp
146	jmc	1.8	CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
147			& SQUEEZE_RIGHT , myThid)
148	jmc	1.9	_END_MASTER( myThid )
149	mlosch	1.1
150	jmc	1.4	C -- end npiter loop
151	jmc	1.3	ENDDO
152	mlosch	1.1	C print some diagnostics
153	jmc	1.9	_BEGIN_MASTER( myThid )
154	jmc	1.10	IF ( rmspp .NE. 0. ) THEN
155	jmc	1.3	IF ( rmspp .EQ. rmsppold ) THEN
156	jmc	1.4	WRITE(msgBuf,'(A)')
157	jmc	1.8	& 'Initial hydrostatic pressure did not converge perfectly,'
158			CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
159			& SQUEEZE_RIGHT , myThid)
160	jmc	1.4	WRITE(msgBuf,'(A)')
161	mlosch	1.1	& 'but the RMS-difference is constant, indicating that the'
162	jmc	1.8	CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
163			& SQUEEZE_RIGHT , myThid)
164	jmc	1.4	WRITE(msgBuf,'(A)')
165	jmc	1.8	& 'algorithm converged within machine precision.'
166			CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
167			& SQUEEZE_RIGHT , myThid)
168	jmc	1.3	ELSE
169	jmc	1.4	WRITE(msgBuf,'(A,I2,A)')
170	jmc	1.8	& 'Initial hydrostatic pressure did not converge after ',
171	mlosch	1.1	& npiter-1, ' steps'
172			CALL PRINT_ERROR( msgBuf, myThid )
173			STOP 'ABNORMAL END: S/R INI_PRESSURE'
174	jmc	1.3	ENDIF
175			ENDIF
176	jmc	1.4	WRITE(msgBuf,'(A)')
177	jmc	1.8	& 'Initial hydrostatic pressure converged.'
178			CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
179			& SQUEEZE_RIGHT , myThid)
180	jmc	1.9	_END_MASTER( myThid )
181	mlosch	1.1
182	heimbach	1.5	#endif /* ALLOW_AUTODIFF_TAMC */
183
184			c-- else of if useDynP_inEos_Zc
185	jmc	1.3	ELSE
186			C print a message and DO nothing
187	jmc	1.9	_BEGIN_MASTER( myThid )
188	jmc	1.4	WRITE(msgBuf,'(A,A)')
189	mlosch	1.1	& 'Pressure is predetermined for buoyancyRelation ',
190			& buoyancyRelation(1:11)
191	jmc	1.8	CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
192			& SQUEEZE_RIGHT , myThid)
193	jmc	1.9	_END_MASTER( myThid )
194	mlosch	1.1
195	jmc	1.3	ENDIF
196	mlosch	1.1
197	jmc	1.9	_BEGIN_MASTER( myThid )
198	jmc	1.4	WRITE(msgBuf,'(A)') ' '
199	jmc	1.8	CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
200			& SQUEEZE_RIGHT , myThid)
201	jmc	1.9	_END_MASTER( myThid )
202	mlosch	1.1
203	jmc	1.8	RETURN
204	jmc	1.3	END