model/src/ini_pressure.F

C $Header: /u/gcmpack/MITgcm/model/src/ini_pressure.F,v 1.15 2013/08/14 18:10:59 jmc Exp $
C $Name:  $

#include "PACKAGES_CONFIG.h"
#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 ( storePhiHyd4Phys ) THEN

#ifndef ALLOW_AUTODIFF
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
        IF ( rmspp.GT.zeroRL ) THEN
         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 )

        ENDIF
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 */

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