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

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

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

revision 1.3 by jmc, Sat Feb 8 02:09:20 2003 UTC revision 1.7 by jmc, Wed Apr 6 18:29:53 2005 UTC
# Line 1  Line 1 
1  C    C $Header$
2  C    C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    
# Line 42  C     bi,bj  - Loop counters Line 42  C     bi,bj  - Loop counters
42  C     I,J,K  C     I,J,K
43        INTEGER bi, bj        INTEGER bi, bj
44        INTEGER  I,  J, K        INTEGER  I,  J, K
45        INTEGER ncount, iMin, iMax, jMin, jMax, npiter        INTEGER  iMin, iMax, jMin, jMax, npiter
46        _RL phiHyd(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL PhiHydF (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
47          _RL PhiHydC (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
48        _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
49        _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
50        _RL pold(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL oldPhi  (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
51        _RL rmspp, rmsppold        _RL count, rmspp, rmsppold
52    
53        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
54  CEOP  CEOP
# Line 60  CEOP Line 61  CEOP
61       &     'Start initial hydrostatic pressure computation'       &     'Start initial hydrostatic pressure computation'
62        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
63       &     SQUEEZE_RIGHT , 1)       &     SQUEEZE_RIGHT , 1)
 C     initialise pressure to water column thickness times a constant  
64        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
65         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
66          DO k = 1,Nr          DO k = 1,Nr
67           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
68            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
69             phiHyd(i,j,k) = 0. _d 0             totPhiHyd(i,j,k,bi,bj) = 0. _d 0
            dPhiHydX(i,j) = 0. _d 0  
            dPhiHydY(i,j) = 0. _d 0  
70            ENDDO            ENDDO
71           ENDDO           ENDDO
          CALL STORE_PRESSURE(  bi, bj, k, phiHyd, myThid )  
72          ENDDO          ENDDO
73         ENDDO         ENDDO
74        ENDDO        ENDDO
75    
76        IF ( buoyancyRelation .EQ. 'OCEANIC' ) THEN        IF ( useDynP_inEos_Zc ) THEN
77    
78    #ifndef ALLOW_AUTODIFF_TAMC
79    cph-- deal with this iterative loop for AD once it will
80    cph-- really be needed;
81    cph-- would need storing of totPhiHyd for each npiter
82    
83         rmspp    = 1. _d 0         rmspp    = 1. _d 0
84         rmsppold = 0. _d 0         rmsppold = 0. _d 0
85         npiter = 0         npiter = 0
86    
87  C     iterate pressure p = integral of (g*rho(p)*dz)  C     iterate pressure p = integral of (g*rho(p)*dz)
88         DO npiter = 1, 15         DO npiter= 1, 15
89          rmsppold = rmspp          rmsppold = rmspp
90          rmspp    = 0. _d 0          rmspp    = 0. _d 0
91          ncount   = 0          count    = 0.
92          DO bj = myByLo(myThid), myByHi(myThid)          DO bj = myByLo(myThid), myByHi(myThid)
93           DO bi = myBxLo(myThid), myBxHi(myThid)           DO bi = myBxLo(myThid), myBxHi(myThid)
94               DO j=1-OLy,sNy+OLy
95                DO i=1-OLx,sNx+OLx
96                 phiHydF(i,j) = 0. _d 0
97                ENDDO
98               ENDDO
99            DO k = 1, Nr            DO k = 1, Nr
100  C     for each level save old pressure and compute new pressure  C     for each level save old pressure and compute new pressure
101             DO j=jMin,jMax             DO j=jMin,jMax
102              DO i=iMin,iMax              DO i=iMin,iMax
103               pold(i,j,k) = pressure(i,j,k,bi,bj)               oldPhi(i,j) = totPhiHyd(i,j,k,bi,bj)
104              ENDDO              ENDDO
105             ENDDO             ENDDO
106             CALL CALC_PHI_HYD(             CALL CALC_PHI_HYD(
107       I          bi, bj, iMin, iMax, jMin, jMax, k,       I          bi, bj, iMin, iMax, jMin, jMax, k,
108       I          theta, salt,       I          theta, salt,
109       U          phiHyd,       U          phiHydF,
110       O          dPhiHydX, dPhiHydY,       O          phiHydC, dPhiHydX, dPhiHydY,
111       I          startTime, nIter0, myThid)       I          startTime, nIter0, myThid)
            CALL STORE_PRESSURE(  bi, bj, k, phiHyd, myThid )  
           ENDDO  
112  C     compute convergence criterion  C     compute convergence criterion
           DO k = 1, Nr  
113             DO j=jMin,jMax             DO j=jMin,jMax
114              DO i=iMin,iMax              DO i=iMin,iMax
115               rmspp = rmspp               rmspp = rmspp
116       &            + (pressure(i,j,k,bi,bj)-pold(i,j,k))**2       &            + (totPhiHyd(i,j,k,bi,bj)-oldPhi(i,j))**2
117               ncount = ncount + int(maskC(i,j,k,bi,bj))               count = count + maskC(i,j,k,bi,bj)
118              ENDDO              ENDDO
119             ENDDO             ENDDO
120    C --      end k loop
121            ENDDO            ENDDO
122           ENDDO           ENDDO
123          ENDDO          ENDDO
124  Cml        WRITE(msgBuf,'(I10.10)') npiter  Cml        WRITE(msgBuf,'(I10.10)') npiter
125  Cml        CALL WRITE_FLD_XYZ_RL( 'POLD.',msgBuf,pold,npiter,myThid)  Cml        CALL WRITE_FLD_XYZ_RL( 'POLD.',msgBuf,pold,npiter,myThid)
126  Cml        CALL WRITE_FLD_XYZ_RL( 'PINI.',msgBuf,pressure,npiter,myThid)  Cml        CALL WRITE_FLD_XYZ_RL( 'PINI.',msgBuf,pressure,npiter,myThid)
127          IF ( ncount .EQ. 0 ) THEN          _GLOBAL_SUM_R8( rmspp, myThid )
128            _GLOBAL_SUM_R8( count, myThid )  
129            IF ( count .EQ. 0. ) THEN
130             rmspp = 0. _d 0             rmspp = 0. _d 0
131          ELSE          ELSE
132             rmspp = sqrt(rmspp/ncount)             rmspp = sqrt(rmspp/count)
133          ENDIF          ENDIF
134          WRITE(msgBuf,'(a,i2,a,e20.13)')          WRITE(msgBuf,'(a,i2,a,e20.13)')
135       &       'Iteration ', npiter, ', RMS-difference = ', rmspp       &       'Iteration ', npiter, ', RMS-difference = ', rmspp
136          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
137       &       SQUEEZE_RIGHT , 1)       &       SQUEEZE_RIGHT , 1)
138    
139    C --   end npiter loop
140         ENDDO         ENDDO
141  C     print some diagnostics      C     print some diagnostics    
142         IF ( rmspp .ne. 0. ) THEN         IF ( rmspp .ne. 0. ) THEN
143          IF ( rmspp .EQ. rmsppold ) THEN          IF ( rmspp .EQ. rmsppold ) THEN
144           WRITE(msgBuf,'(a)')           WRITE(msgBuf,'(A)')
145       &      'Initial hydrostatic pressure did not converge perfectly,'       &      'Initial hydrostatic pressure did not converge perfectly,'
146           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
147       &        SQUEEZE_RIGHT , 1)       &        SQUEEZE_RIGHT , 1)
148           WRITE(msgBuf,'(a)')           WRITE(msgBuf,'(A)')
149       &      'but the RMS-difference is constant, indicating that the'       &      'but the RMS-difference is constant, indicating that the'
150           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
151       &        SQUEEZE_RIGHT , 1)       &        SQUEEZE_RIGHT , 1)
152           WRITE(msgBuf,'(a)')           WRITE(msgBuf,'(A)')
153       &      'algorithm converged within machine precision.'           &      'algorithm converged within machine precision.'    
154           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
155       &        SQUEEZE_RIGHT , 1)       &        SQUEEZE_RIGHT , 1)
156          ELSE          ELSE
157           WRITE(msgBuf,'(a,i2,a)')           WRITE(msgBuf,'(A,I2,A)')
158       &       'Initial hydrostatic pressure did not converge after ',       &       'Initial hydrostatic pressure did not converge after ',
159       &          npiter-1, ' steps'       &          npiter-1, ' steps'
160           CALL PRINT_ERROR( msgBuf, myThid )           CALL PRINT_ERROR( msgBuf, myThid )
161           STOP 'ABNORMAL END: S/R INI_PRESSURE'           STOP 'ABNORMAL END: S/R INI_PRESSURE'
162          ENDIF          ENDIF
163         ENDIF         ENDIF
164         WRITE(msgBuf,'(a)')         WRITE(msgBuf,'(A)')
165       &     'Initial hydrostatic pressure converged.'       &     'Initial hydrostatic pressure converged.'
166         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
167       &      SQUEEZE_RIGHT , 1)       &      SQUEEZE_RIGHT , 1)
168    
169    #endif /* ALLOW_AUTODIFF_TAMC */
170    
171    c-- else of if useDynP_inEos_Zc
172        ELSE        ELSE
173  C     print a message and DO nothing  C     print a message and DO nothing
174         WRITE(msgBuf,'(a,a)')         WRITE(msgBuf,'(A,A)')
175       &        'Pressure is predetermined for buoyancyRelation ',       &        'Pressure is predetermined for buoyancyRelation ',
176       &        buoyancyRelation(1:11)       &        buoyancyRelation(1:11)
177         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
# Line 168  C     print a message and DO nothing Line 179  C     print a message and DO nothing
179    
180        ENDIF        ENDIF
181    
182        WRITE(msgBuf,'(a)') ' '        WRITE(msgBuf,'(A)') ' '
183        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
184       &     SQUEEZE_RIGHT , 1)       &     SQUEEZE_RIGHT , 1)
185    
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.22