/[MITgcm]/MITgcm/model/src/ini_pressure.F
ViewVC logotype

Annotation of /MITgcm/model/src/ini_pressure.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.14 - (hide annotations) (download)
Sun Aug 12 20:25:33 2012 UTC (11 years, 9 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64, checkpoint63r, checkpoint63s, checkpoint64i, checkpoint64h, checkpoint64k, checkpoint64j, checkpoint64m, checkpoint64l, checkpoint64a, checkpoint64c, checkpoint64b, checkpoint64e, checkpoint64d, checkpoint64g, checkpoint64f
Changes since 1.13: +4 -3 lines
- add PACKAGES_CONFIG.h wherever ALLOW_AUTODIFF[_TAMC] is used (in model/src);
- also change #ifdef ALLOW_AUTODIFF_TAMC to #ifdef ALLOW_AUTODIFF

1 jmc 1.14 C $Header: /u/gcmpack/MITgcm/model/src/ini_pressure.F,v 1.13 2010/08/12 21:49:55 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     _RL PhiHydF (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
48     _RL PhiHydC (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
49 jmc 1.3 _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
50     _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
51 jmc 1.4 _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     & 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.4 IF ( useDynP_inEos_Zc ) 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 (g*rho(p)*dz)
93 jmc 1.4 DO npiter= 1, 15
94 mlosch 1.1 rmsppold = rmspp
95     rmspp = 0. _d 0
96 jmc 1.10 count = 0. _d 0
97 jmc 1.3 DO bj = myByLo(myThid), myByHi(myThid)
98     DO bi = myBxLo(myThid), myBxHi(myThid)
99 jmc 1.10 rmsTile = 0. _d 0
100     sumTile = 0. _d 0
101     DO j=1-OLy,sNy+OLy
102     DO i=1-OLx,sNx+OLx
103 jmc 1.4 phiHydF(i,j) = 0. _d 0
104     ENDDO
105 jmc 1.10 ENDDO
106 jmc 1.3 DO k = 1, Nr
107 jmc 1.10 C for each level save old pressure and compute new pressure
108 jmc 1.3 DO j=jMin,jMax
109     DO i=iMin,iMax
110 jmc 1.4 oldPhi(i,j) = totPhiHyd(i,j,k,bi,bj)
111 jmc 1.3 ENDDO
112     ENDDO
113     CALL CALC_PHI_HYD(
114     I bi, bj, iMin, iMax, jMin, jMax, k,
115 jmc 1.8 I theta, salt,
116 jmc 1.4 U phiHydF,
117     O phiHydC, dPhiHydX, dPhiHydY,
118 jmc 1.11 I startTime, -1, myThid )
119 mlosch 1.1 C compute convergence criterion
120 jmc 1.10 DO j=1,sNy
121     DO i=1,sNx
122     rmsTile = rmsTile
123 jmc 1.4 & + (totPhiHyd(i,j,k,bi,bj)-oldPhi(i,j))**2
124 jmc 1.10 & * maskC(i,j,k,bi,bj)
125     sumTile = sumTile + maskC(i,j,k,bi,bj)
126 jmc 1.3 ENDDO
127     ENDDO
128 jmc 1.4 C -- end k loop
129 jmc 1.3 ENDDO
130 jmc 1.10 rmspp = rmspp + rmsTile
131     count = count + sumTile
132 jmc 1.3 ENDDO
133     ENDDO
134 mlosch 1.1 Cml WRITE(msgBuf,'(I10.10)') npiter
135     Cml CALL WRITE_FLD_XYZ_RL( 'POLD.',msgBuf,pold,npiter,myThid)
136     Cml CALL WRITE_FLD_XYZ_RL( 'PINI.',msgBuf,pressure,npiter,myThid)
137 jmc 1.12 _GLOBAL_SUM_RL( rmspp, myThid )
138     _GLOBAL_SUM_RL( count, myThid )
139 jmc 1.4 IF ( count .EQ. 0. ) THEN
140 mlosch 1.1 rmspp = 0. _d 0
141 jmc 1.3 ELSE
142 jmc 1.10 rmspp = SQRT(rmspp/count)
143 jmc 1.3 ENDIF
144 jmc 1.9 _BEGIN_MASTER( myThid )
145 jmc 1.10 WRITE(msgBuf,'(A,I4,A,1PE20.12)')
146     & 'Iteration', npiter, ', RMS-difference =', rmspp
147 jmc 1.8 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
148     & SQUEEZE_RIGHT , myThid)
149 jmc 1.9 _END_MASTER( myThid )
150 mlosch 1.1
151 jmc 1.4 C -- end npiter loop
152 jmc 1.3 ENDDO
153 jmc 1.13 C print some diagnostics
154 jmc 1.9 _BEGIN_MASTER( myThid )
155 jmc 1.10 IF ( rmspp .NE. 0. ) THEN
156 jmc 1.3 IF ( rmspp .EQ. rmsppold ) THEN
157 jmc 1.4 WRITE(msgBuf,'(A)')
158 jmc 1.8 & 'Initial hydrostatic pressure did not converge perfectly,'
159     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
160     & SQUEEZE_RIGHT , myThid)
161 jmc 1.4 WRITE(msgBuf,'(A)')
162 mlosch 1.1 & 'but the RMS-difference is constant, indicating that the'
163 jmc 1.8 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
164     & SQUEEZE_RIGHT , myThid)
165 jmc 1.4 WRITE(msgBuf,'(A)')
166 jmc 1.8 & 'algorithm converged within machine precision.'
167     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
168     & SQUEEZE_RIGHT , myThid)
169 jmc 1.3 ELSE
170 jmc 1.4 WRITE(msgBuf,'(A,I2,A)')
171 jmc 1.8 & 'Initial hydrostatic pressure did not converge after ',
172 mlosch 1.1 & npiter-1, ' steps'
173     CALL PRINT_ERROR( msgBuf, myThid )
174     STOP 'ABNORMAL END: S/R INI_PRESSURE'
175 jmc 1.3 ENDIF
176     ENDIF
177 jmc 1.4 WRITE(msgBuf,'(A)')
178 jmc 1.8 & 'Initial hydrostatic pressure converged.'
179     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
180     & SQUEEZE_RIGHT , myThid)
181 jmc 1.9 _END_MASTER( myThid )
182 mlosch 1.1
183 jmc 1.14 #endif /* ALLOW_AUTODIFF */
184 heimbach 1.5
185     c-- else of if useDynP_inEos_Zc
186 jmc 1.3 ELSE
187     C print a message and DO nothing
188 jmc 1.9 _BEGIN_MASTER( myThid )
189 jmc 1.4 WRITE(msgBuf,'(A,A)')
190 mlosch 1.1 & 'Pressure is predetermined for buoyancyRelation ',
191     & buoyancyRelation(1:11)
192 jmc 1.8 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
193     & SQUEEZE_RIGHT , myThid)
194 jmc 1.9 _END_MASTER( myThid )
195 mlosch 1.1
196 jmc 1.3 ENDIF
197 mlosch 1.1
198 jmc 1.9 _BEGIN_MASTER( myThid )
199 jmc 1.4 WRITE(msgBuf,'(A)') ' '
200 jmc 1.8 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
201     & SQUEEZE_RIGHT , myThid)
202 jmc 1.9 _END_MASTER( myThid )
203 mlosch 1.1
204 jmc 1.8 RETURN
205 jmc 1.3 END

  ViewVC Help
Powered by ViewVC 1.1.22