/[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.2 by mlosch, Wed Sep 25 19:36:50 2002 UTC revision 1.3 by jmc, Sat Feb 8 02:09:20 2003 UTC
# Line 6  C Line 6  C
6  CBOP  CBOP
7  C     !ROUTINE: INI_PRESSURE  C     !ROUTINE: INI_PRESSURE
8  C     !INTERFACE:  C     !INTERFACE:
9        subroutine ini_pressure( myThid )        SUBROUTINE INI_PRESSURE( myThid )
10  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
11  C     *==========================================================*  C     *==========================================================*
12  C     | SUBROUTINE INI_PRESSURE  C     | SUBROUTINE INI_PRESSURE
# Line 16  C     |   - needs to be called after ini Line 16  C     |   - needs to be called after ini
16  C     |     ini_psurf  C     |     ini_psurf
17  C     |   - does not include surface pressure loading, because  C     |   - does not include surface pressure loading, because
18  C     |     that is only available until after  C     |     that is only available until after
19  C     |     call packages_init_variables  C     |     CALL packages_init_variables
20  C     *==========================================================*  C     *==========================================================*
21  C     \ev  C     \ev
22    
23  C     !USES:  C     !USES:
24        implicit none        IMPLICIT NONE
25  C     == Global variables ==  C     == Global variables ==
26  #include "SIZE.h"  #include "SIZE.h"
27  #include "EEPARAMS.h"  #include "EEPARAMS.h"
# Line 37  C     myThid -  Number of this instance Line 37  C     myThid -  Number of this instance
37    
38  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
39  C     == Local variables ==  C     == Local variables ==
40    C     dPhiHydX,Y :: Gradient (X & Y directions) of Hyd. Potential
41  C     bi,bj  - Loop counters  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 ncount, iMin, iMax, jMin, jMax, npiter
46        _RL phiHyd(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL phiHyd(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
47          _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
48          _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
49        _RL pold(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL pold(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
50        _RL rmspp, rmsppold        _RL rmspp, rmsppold
51    
# Line 53  CEOP Line 56  CEOP
56        jMin = 1-OLy        jMin = 1-OLy
57        jMax = sNy+OLy        jMax = sNy+OLy
58                
59        write(msgBuf,'(a)')        WRITE(msgBuf,'(a)')
60       &     'Start initial hydrostatic pressure computation'       &     'Start initial hydrostatic pressure computation'
61        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
62       &     SQUEEZE_RIGHT , 1)       &     SQUEEZE_RIGHT , 1)
63  C     initialise pressure to water column thickness times a constant  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             phiHyd(i,j,k) = 0. _d 0
70            end do             dPhiHydX(i,j) = 0. _d 0
71           end do             dPhiHydY(i,j) = 0. _d 0
72           call store_pressure(  bi, bj, k, phiHyd, myThid )            ENDDO
73          end do           ENDDO
74         end do           CALL STORE_PRESSURE(  bi, bj, k, phiHyd, myThid )
75        end do          ENDDO
76           ENDDO
77          ENDDO
78    
79        if ( buoyancyRelation .eq. 'OCEANIC' ) then        IF ( buoyancyRelation .EQ. 'OCEANIC' ) THEN
80    
81         rmspp    = 1. _d 0         rmspp    = 1. _d 0
82         rmsppold = 0. _d 0         rmsppold = 0. _d 0
83         npiter = 0         npiter = 0
84    
85  C     iterate pressure p = integral of (g*rho(p)*dz)  C     iterate pressure p = integral of (g*rho(p)*dz)
86         do npiter = 1, 15         DO npiter = 1, 15
87          rmsppold = rmspp          rmsppold = rmspp
88          rmspp    = 0. _d 0          rmspp    = 0. _d 0
89          ncount   = 0          ncount   = 0
90          do bj = myByLo(myThid), myByHi(myThid)          DO bj = myByLo(myThid), myByHi(myThid)
91           do bi = myBxLo(myThid), myBxHi(myThid)           DO bi = myBxLo(myThid), myBxHi(myThid)
92            do k = 1, Nr            DO k = 1, Nr
93  C     for each level save old pressure and compute new pressure  C     for each level save old pressure and compute new pressure
94             do j=jMin,jMax             DO j=jMin,jMax
95              do i=iMin,iMax              DO i=iMin,iMax
96               pold(i,j,k) = pressure(i,j,k,bi,bj)               pold(i,j,k) = pressure(i,j,k,bi,bj)
97              end do              ENDDO
98             end do             ENDDO
99             call calc_phi_hyd(             CALL CALC_PHI_HYD(
100       &          bi, bj, iMin, iMax, jMin, jMax, k,       I          bi, bj, iMin, iMax, jMin, jMax, k,
101       &          theta, salt, phiHyd, myThid)       I          theta, salt,
102             call store_pressure(  bi, bj, k, phiHyd, myThid )       U          phiHyd,
103            end do       O          dPhiHydX, dPhiHydY,
104         I          startTime, nIter0, myThid)
105               CALL STORE_PRESSURE(  bi, bj, k, phiHyd, myThid )
106              ENDDO
107  C     compute convergence criterion  C     compute convergence criterion
108            do k = 1, Nr            DO k = 1, Nr
109             do j=jMin,jMax             DO j=jMin,jMax
110              do i=iMin,iMax              DO i=iMin,iMax
111               rmspp = rmspp               rmspp = rmspp
112       &            + (pressure(i,j,k,bi,bj)-pold(i,j,k))**2       &            + (pressure(i,j,k,bi,bj)-pold(i,j,k))**2
113               ncount = ncount + int(maskC(i,j,k,bi,bj))               ncount = ncount + int(maskC(i,j,k,bi,bj))
114              end do              ENDDO
115             end do             ENDDO
116            end do            ENDDO
117           end do           ENDDO
118          end do          ENDDO
119  Cml        WRITE(msgBuf,'(I10.10)') npiter  Cml        WRITE(msgBuf,'(I10.10)') npiter
120  Cml        CALL WRITE_FLD_XYZ_RL( 'POLD.',msgBuf,pold,npiter,myThid)  Cml        CALL WRITE_FLD_XYZ_RL( 'POLD.',msgBuf,pold,npiter,myThid)
121  Cml        CALL WRITE_FLD_XYZ_RL( 'PINI.',msgBuf,pressure,npiter,myThid)  Cml        CALL WRITE_FLD_XYZ_RL( 'PINI.',msgBuf,pressure,npiter,myThid)
122          if ( ncount .eq. 0 ) then          IF ( ncount .EQ. 0 ) THEN
123             rmspp = 0. _d 0             rmspp = 0. _d 0
124          else          ELSE
125             rmspp = sqrt(rmspp/ncount)             rmspp = sqrt(rmspp/ncount)
126          end if          ENDIF
127          write(msgBuf,'(a,i2,a,e20.13)')          WRITE(msgBuf,'(a,i2,a,e20.13)')
128       &       'Iteration ', npiter, ', RMS-difference = ', rmspp       &       'Iteration ', npiter, ', RMS-difference = ', rmspp
129          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
130       &       SQUEEZE_RIGHT , 1)       &       SQUEEZE_RIGHT , 1)
131    
132         end do         ENDDO
133  C     print some diagnostics      C     print some diagnostics    
134         if ( rmspp .ne. 0. ) then         IF ( rmspp .ne. 0. ) THEN
135          if ( rmspp .eq. rmsppold ) then          IF ( rmspp .EQ. rmsppold ) THEN
136           write(msgBuf,'(a)')           WRITE(msgBuf,'(a)')
137       &      'Initial hydrostatic pressure did not converge perfectly,'       &      'Initial hydrostatic pressure did not converge perfectly,'
138           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
139       &        SQUEEZE_RIGHT , 1)       &        SQUEEZE_RIGHT , 1)
140           write(msgBuf,'(a)')           WRITE(msgBuf,'(a)')
141       &      'but the RMS-difference is constant, indicating that the'       &      'but the RMS-difference is constant, indicating that the'
142           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
143       &        SQUEEZE_RIGHT , 1)       &        SQUEEZE_RIGHT , 1)
144           write(msgBuf,'(a)')           WRITE(msgBuf,'(a)')
145       &      'algorithm converged within machine precision.'           &      'algorithm converged within machine precision.'    
146           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
147       &        SQUEEZE_RIGHT , 1)       &        SQUEEZE_RIGHT , 1)
148          else          ELSE
149           write(msgBuf,'(a,i2,a)')           WRITE(msgBuf,'(a,i2,a)')
150       &       'Initial hydrostatic pressure did not converge after ',       &       'Initial hydrostatic pressure did not converge after ',
151       &          npiter-1, ' steps'       &          npiter-1, ' steps'
152           CALL PRINT_ERROR( msgBuf, myThid )           CALL PRINT_ERROR( msgBuf, myThid )
153           STOP 'ABNORMAL END: S/R INI_PRESSURE'           STOP 'ABNORMAL END: S/R INI_PRESSURE'
154          end if          ENDIF
155         end if         ENDIF
156         write(msgBuf,'(a)')         WRITE(msgBuf,'(a)')
157       &     'Initial hydrostatic pressure converged.'       &     'Initial hydrostatic pressure converged.'
158         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
159       &      SQUEEZE_RIGHT , 1)       &      SQUEEZE_RIGHT , 1)
160    
161        else        ELSE
162  C     print a message and do nothing  C     print a message and DO nothing
163         write(msgBuf,'(a,a)')         WRITE(msgBuf,'(a,a)')
164       &        'Pressure is predetermined for buoyancyRelation ',       &        'Pressure is predetermined for buoyancyRelation ',
165       &        buoyancyRelation(1:11)       &        buoyancyRelation(1:11)
166         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
167       &      SQUEEZE_RIGHT , 1)       &      SQUEEZE_RIGHT , 1)
168    
169        end if        ENDIF
170    
171        write(msgBuf,'(a)') ' '        WRITE(msgBuf,'(a)') ' '
172        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
173       &     SQUEEZE_RIGHT , 1)       &     SQUEEZE_RIGHT , 1)
174    
175        return        RETURN
176        end        END
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.22