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

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

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

Revision 1.2 - (show 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 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 CEOP
51 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