/[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.4 by jmc, Tue Feb 18 15:34:25 2003 UTC
# 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  c     IF ( startTime .NE. 0. ) RETURN
77    c     IF ( buoyancyRelation .EQ. 'OCEANIC' ) THEN
78          IF ( useDynP_inEos_Zc ) THEN
79    
80         rmspp    = 1. _d 0         rmspp    = 1. _d 0
81         rmsppold = 0. _d 0         rmsppold = 0. _d 0
82         npiter = 0         npiter = 0
83    
84  C     iterate pressure p = integral of (g*rho(p)*dz)  C     iterate pressure p = integral of (g*rho(p)*dz)
85         DO npiter = 1, 15         DO npiter= 1, 15
86          rmsppold = rmspp          rmsppold = rmspp
87          rmspp    = 0. _d 0          rmspp    = 0. _d 0
88          ncount   = 0          count    = 0.
89          DO bj = myByLo(myThid), myByHi(myThid)          DO bj = myByLo(myThid), myByHi(myThid)
90           DO bi = myBxLo(myThid), myBxHi(myThid)           DO bi = myBxLo(myThid), myBxHi(myThid)
91               DO j=1-OLy,sNy+OLy
92                DO i=1-OLx,sNx+OLx
93                 phiHydF(i,j) = 0. _d 0
94                ENDDO
95               ENDDO
96            DO k = 1, Nr            DO k = 1, Nr
97  C     for each level save old pressure and compute new pressure  C     for each level save old pressure and compute new pressure
98             DO j=jMin,jMax             DO j=jMin,jMax
99              DO i=iMin,iMax              DO i=iMin,iMax
100               pold(i,j,k) = pressure(i,j,k,bi,bj)               oldPhi(i,j) = totPhiHyd(i,j,k,bi,bj)
101              ENDDO              ENDDO
102             ENDDO             ENDDO
103             CALL CALC_PHI_HYD(             CALL CALC_PHI_HYD(
104       I          bi, bj, iMin, iMax, jMin, jMax, k,       I          bi, bj, iMin, iMax, jMin, jMax, k,
105       I          theta, salt,       I          theta, salt,
106       U          phiHyd,       U          phiHydF,
107       O          dPhiHydX, dPhiHydY,       O          phiHydC, dPhiHydX, dPhiHydY,
108       I          startTime, nIter0, myThid)       I          startTime, nIter0, myThid)
            CALL STORE_PRESSURE(  bi, bj, k, phiHyd, myThid )  
           ENDDO  
109  C     compute convergence criterion  C     compute convergence criterion
           DO k = 1, Nr  
110             DO j=jMin,jMax             DO j=jMin,jMax
111              DO i=iMin,iMax              DO i=iMin,iMax
112               rmspp = rmspp               rmspp = rmspp
113       &            + (pressure(i,j,k,bi,bj)-pold(i,j,k))**2       &            + (totPhiHyd(i,j,k,bi,bj)-oldPhi(i,j))**2
114               ncount = ncount + int(maskC(i,j,k,bi,bj))               count = count + maskC(i,j,k,bi,bj)
115              ENDDO              ENDDO
116             ENDDO             ENDDO
117    C --      end k loop
118            ENDDO            ENDDO
119           ENDDO           ENDDO
120          ENDDO          ENDDO
121  Cml        WRITE(msgBuf,'(I10.10)') npiter  Cml        WRITE(msgBuf,'(I10.10)') npiter
122  Cml        CALL WRITE_FLD_XYZ_RL( 'POLD.',msgBuf,pold,npiter,myThid)  Cml        CALL WRITE_FLD_XYZ_RL( 'POLD.',msgBuf,pold,npiter,myThid)
123  Cml        CALL WRITE_FLD_XYZ_RL( 'PINI.',msgBuf,pressure,npiter,myThid)  Cml        CALL WRITE_FLD_XYZ_RL( 'PINI.',msgBuf,pressure,npiter,myThid)
124          IF ( ncount .EQ. 0 ) THEN          _GLOBAL_SUM_R8( rmspp, myThid )
125            _GLOBAL_SUM_R8( count, myThid )  
126            IF ( count .EQ. 0. ) THEN
127             rmspp = 0. _d 0             rmspp = 0. _d 0
128          ELSE          ELSE
129             rmspp = sqrt(rmspp/ncount)             rmspp = sqrt(rmspp/count)
130          ENDIF          ENDIF
131          WRITE(msgBuf,'(a,i2,a,e20.13)')          WRITE(msgBuf,'(a,i2,a,e20.13)')
132       &       'Iteration ', npiter, ', RMS-difference = ', rmspp       &       'Iteration ', npiter, ', RMS-difference = ', rmspp
133          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
134       &       SQUEEZE_RIGHT , 1)       &       SQUEEZE_RIGHT , 1)
135    
136    C --   end npiter loop
137         ENDDO         ENDDO
138  C     print some diagnostics      C     print some diagnostics    
139         IF ( rmspp .ne. 0. ) THEN         IF ( rmspp .ne. 0. ) THEN
140          IF ( rmspp .EQ. rmsppold ) THEN          IF ( rmspp .EQ. rmsppold ) THEN
141           WRITE(msgBuf,'(a)')           WRITE(msgBuf,'(A)')
142       &      'Initial hydrostatic pressure did not converge perfectly,'       &      'Initial hydrostatic pressure did not converge perfectly,'
143           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
144       &        SQUEEZE_RIGHT , 1)       &        SQUEEZE_RIGHT , 1)
145           WRITE(msgBuf,'(a)')           WRITE(msgBuf,'(A)')
146       &      'but the RMS-difference is constant, indicating that the'       &      'but the RMS-difference is constant, indicating that the'
147           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
148       &        SQUEEZE_RIGHT , 1)       &        SQUEEZE_RIGHT , 1)
149           WRITE(msgBuf,'(a)')           WRITE(msgBuf,'(A)')
150       &      'algorithm converged within machine precision.'           &      'algorithm converged within machine precision.'    
151           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
152       &        SQUEEZE_RIGHT , 1)       &        SQUEEZE_RIGHT , 1)
153          ELSE          ELSE
154           WRITE(msgBuf,'(a,i2,a)')           WRITE(msgBuf,'(A,I2,A)')
155       &       'Initial hydrostatic pressure did not converge after ',       &       'Initial hydrostatic pressure did not converge after ',
156       &          npiter-1, ' steps'       &          npiter-1, ' steps'
157           CALL PRINT_ERROR( msgBuf, myThid )           CALL PRINT_ERROR( msgBuf, myThid )
158           STOP 'ABNORMAL END: S/R INI_PRESSURE'           STOP 'ABNORMAL END: S/R INI_PRESSURE'
159          ENDIF          ENDIF
160         ENDIF         ENDIF
161         WRITE(msgBuf,'(a)')         WRITE(msgBuf,'(A)')
162       &     'Initial hydrostatic pressure converged.'       &     'Initial hydrostatic pressure converged.'
163         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
164       &      SQUEEZE_RIGHT , 1)       &      SQUEEZE_RIGHT , 1)
165    
166        ELSE        ELSE
167  C     print a message and DO nothing  C     print a message and DO nothing
168         WRITE(msgBuf,'(a,a)')         WRITE(msgBuf,'(A,A)')
169       &        'Pressure is predetermined for buoyancyRelation ',       &        'Pressure is predetermined for buoyancyRelation ',
170       &        buoyancyRelation(1:11)       &        buoyancyRelation(1:11)
171         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
# Line 168  C     print a message and DO nothing Line 173  C     print a message and DO nothing
173    
174        ENDIF        ENDIF
175    
176        WRITE(msgBuf,'(a)') ' '        WRITE(msgBuf,'(A)') ' '
177        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
178       &     SQUEEZE_RIGHT , 1)       &     SQUEEZE_RIGHT , 1)
179    
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.22