/[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.2 - (hide annotations) (download)
Wed Sep 25 19:36:50 2002 UTC (21 years, 8 months ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint46n_post, checkpoint47e_post, checkpoint46l_post, checkpoint47c_post, checkpoint46l_pre, checkpoint48b_post, checkpoint48c_pre, checkpoint47d_pre, checkpoint47a_post, checkpoint47i_post, checkpoint47d_post, checkpoint46j_pre, checkpoint47g_post, checkpoint46j_post, checkpoint46k_post, checkpoint48a_post, checkpoint47j_post, branch-exfmods-tag, checkpoint48c_post, checkpoint47b_post, checkpoint46m_post, checkpoint47f_post, checkpoint46i_post, checkpoint47, checkpoint48, checkpoint46h_post, checkpoint47h_post
Branch point for: branch-exfmods-curt
Changes since 1.1: +1 -1 lines 
o cleaned up the use of rhoNil and rhoConst.
  - rhoNil should only appear in the LINEAR equation of state, everywhere
    else rhoNil is replaced by rhoConst, e.g. find_rho computes rho-rhoConst
    and the dynamical equations are all divided by rhoConst
o introduced new parameter rhoConstFresh, a reference density of fresh
  water, to remove the fresh water flux's dependence on rhoNil. The default
  value is 999.8 kg/m^3
o cleanup up external_forcing.F and external_forcing_surf.F
  - can now be used by both OCEANIC and OCEANICP
1 mlosch 1.1 C $Header:
2     C $Name:
3    
4     #include "CPP_OPTIONS.h"
5    
6     CBOP
7     C !ROUTINE: INI_PRESSURE
8     C !INTERFACE:
9     subroutine ini_pressure( myThid )
10     C !DESCRIPTION: \bv
11     C *==========================================================*
12     C | SUBROUTINE INI_PRESSURE
13     C | o initialise the pressure field consistently with
14     C | temperature and salinity
15     C | - needs to be called after ini_theta, ini_salt, and
16     C | ini_psurf
17     C | - does not include surface pressure loading, because
18     C | that is only available until after
19     C | call packages_init_variables
20     C *==========================================================*
21     C \ev
22    
23     C !USES:
24     implicit none
25     C == Global variables ==
26     #include "SIZE.h"
27     #include "EEPARAMS.h"
28     #include "PARAMS.h"
29     #include "EOS.h"
30     #include "GRID.h"
31     #include "DYNVARS.h"
32    
33     C !INPUT/OUTPUT PARAMETERS:
34     C == Routine arguments ==
35     C myThid - Number of this instance of INI_PRESSURE
36     INTEGER myThid
37    
38     C !LOCAL VARIABLES:
39     C == Local variables ==
40     C bi,bj - Loop counters
41     C I,J,K
42     INTEGER bi, bj
43     INTEGER I, J, K
44     INTEGER ncount, iMin, iMax, jMin, jMax, npiter
45     _RL phiHyd(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
46     _RL pold(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
47     _RL rmspp, rmsppold
48    
49     CHARACTER*(MAX_LEN_MBUF) msgBuf
50 mlosch 1.2 CEOP
51 mlosch 1.1 iMin = 1-OLx
52     iMax = sNx+OLx
53     jMin = 1-OLy
54     jMax = sNy+OLy
55    
56     write(msgBuf,'(a)')
57     & 'Start initial hydrostatic pressure computation'
58     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
59     & SQUEEZE_RIGHT , 1)
60     C initialise pressure to water column thickness times a constant
61     do bj = myByLo(myThid), myByHi(myThid)
62     do bi = myBxLo(myThid), myBxHi(myThid)
63     do k = 1,Nr
64     do j=1-OLy,sNy+OLy
65     do i=1-OLx,sNx+OLx
66     phiHyd(i,j,k) = 0. _d 0
67     end do
68     end do
69     call store_pressure( bi, bj, k, phiHyd, myThid )
70     end do
71     end do
72     end do
73    
74     if ( buoyancyRelation .eq. 'OCEANIC' ) then
75    
76     rmspp = 1. _d 0
77     rmsppold = 0. _d 0
78     npiter = 0
79    
80     C iterate pressure p = integral of (g*rho(p)*dz)
81     do npiter = 1, 15
82     rmsppold = rmspp
83     rmspp = 0. _d 0
84     ncount = 0
85     do bj = myByLo(myThid), myByHi(myThid)
86     do bi = myBxLo(myThid), myBxHi(myThid)
87     do k = 1, Nr
88     C for each level save old pressure and compute new pressure
89     do j=jMin,jMax
90     do i=iMin,iMax
91     pold(i,j,k) = pressure(i,j,k,bi,bj)
92     end do
93     end do
94     call calc_phi_hyd(
95     & bi, bj, iMin, iMax, jMin, jMax, k,
96     & theta, salt, phiHyd, myThid)
97     call store_pressure( bi, bj, k, phiHyd, myThid )
98     end do
99     C compute convergence criterion
100     do k = 1, Nr
101     do j=jMin,jMax
102     do i=iMin,iMax
103     rmspp = rmspp
104     & + (pressure(i,j,k,bi,bj)-pold(i,j,k))**2
105     ncount = ncount + int(maskC(i,j,k,bi,bj))
106     end do
107     end do
108     end do
109     end do
110     end do
111     Cml WRITE(msgBuf,'(I10.10)') npiter
112     Cml CALL WRITE_FLD_XYZ_RL( 'POLD.',msgBuf,pold,npiter,myThid)
113     Cml CALL WRITE_FLD_XYZ_RL( 'PINI.',msgBuf,pressure,npiter,myThid)
114     if ( ncount .eq. 0 ) then
115     rmspp = 0. _d 0
116     else
117     rmspp = sqrt(rmspp/ncount)
118     end if
119     write(msgBuf,'(a,i2,a,e20.13)')
120     & 'Iteration ', npiter, ', RMS-difference = ', rmspp
121     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
122     & SQUEEZE_RIGHT , 1)
123    
124     end do
125     C print some diagnostics
126     if ( rmspp .ne. 0. ) then
127     if ( rmspp .eq. rmsppold ) then
128     write(msgBuf,'(a)')
129     & 'Initial hydrostatic pressure did not converge perfectly,'
130     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
131     & SQUEEZE_RIGHT , 1)
132     write(msgBuf,'(a)')
133     & 'but the RMS-difference is constant, indicating that the'
134     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
135     & SQUEEZE_RIGHT , 1)
136     write(msgBuf,'(a)')
137     & 'algorithm converged within machine precision.'
138     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
139     & SQUEEZE_RIGHT , 1)
140     else
141     write(msgBuf,'(a,i2,a)')
142     & 'Initial hydrostatic pressure did not converge after ',
143     & npiter-1, ' steps'
144     CALL PRINT_ERROR( msgBuf, myThid )
145     STOP 'ABNORMAL END: S/R INI_PRESSURE'
146     end if
147     end if
148     write(msgBuf,'(a)')
149     & 'Initial hydrostatic pressure converged.'
150     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
151     & SQUEEZE_RIGHT , 1)
152    
153     else
154     C print a message and do nothing
155     write(msgBuf,'(a,a)')
156     & 'Pressure is predetermined for buoyancyRelation ',
157     & buoyancyRelation(1:11)
158     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
159     & SQUEEZE_RIGHT , 1)
160    
161     end if
162    
163     write(msgBuf,'(a)') ' '
164     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
165     & SQUEEZE_RIGHT , 1)
166    
167     return
168     end
  ViewVC Help
Powered by ViewVC 1.1.22