intpressure.F90 5.63 KB
Newer Older
dumoda01's avatar
dumoda01 committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
!$Id: intpressure.F90,v 1.8 2007-01-06 11:49:15 kbk Exp $
#include"cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: The internal pressure-gradient \label{sec:intpressure}
!
! !INTERFACE:
   subroutine intpressure(nlev)
!
! !DESCRIPTION:
!   With the hydrostatic assumption
!  \begin{equation}\label{HydroStat}
!   \partder{P}{z} + g \mean{\rho} = 0
!   \comma
!  \end{equation}
!  where $P$ denotes the mean pressure, $g=9.81$m\,s$^{-2}$
!  the gravitational acceleration  and $\mean{\rho}$ the mean density,
!  the components of the pressure-gradient may be expressed as
!  \begin{equation}
!   \label{InternalPressurex}
!  - \frac{1}{\rho_0} \partder{P}{x}=
!  -g \partder{\zeta}{x}
!  +\int_z^{\zeta}\partder{B}{x} \, dz'
!  \end{equation}
!  and
!  \begin{equation}\label{InternalPressurey}
!  - \frac{1}{\rho_0} \partder{P}{y}=
!  -g \partder{\zeta}{y}
!  +\int_z^{\zeta} \partder{B}{y} \, dz'
!   \comma
!  \end{equation}
!  where $\zeta$ is the surface elevation and $B$ the
!  mean buoyancy as defined in \eq{DefBuoyancy}.
!
!  The first term on the right hand side
!  in \eq{InternalPressurex}
!  and \eq{InternalPressurey} is the external pressure-gradient
!  due to surface slopes,  and the second the internal pressure-gradient
!  due to the density gradient.
!  The internal pressure-gradient will only be established by
!  gradients of mean potential temperature $\Theta$ and mean
!  salinity $S$. Sediment concentration is assumed to be
!  horizontally homogeneous.
!
!  In this subroutine, first, the horizontal buoyancy gradients,
!  $\partial_xB$ and $\partial_yB$,
!  are calculated from the prescribed gradients of salinity, $\partial_xS$
!  and $\partial_yS$, and temperature, $\partial_x\Theta$ and $\partial_y\Theta$,
!  according to the finite-difference expression
!  \begin{equation}
!    \partder{B}{x} \approx
!    \frac{B(S+\Delta_xS,\Theta+\Delta_x\Theta,P)-B(S,\Theta,P)}{\Delta x}
!    \comma
!  \end{equation}
!  \begin{equation}
!    \partder{B}{y} \approx
!    \frac{B(S+\Delta_yS,\Theta+\Delta_y\Theta,P)-B(S,\theta,P)}{\Delta y}
!   \comma
!  \end{equation}
!  where the defintions
!  \begin{equation}
!    \Delta_xS=\Delta x \partial_xS \comma
!    \Delta_yS=\Delta y \partial_yS \comma
!  \end{equation}
!  and
!  \begin{equation}
!   \Delta_x\Theta=\Delta x \partial_x\Theta \comma
!   \Delta_y\Theta=\Delta y \partial_y\Theta \comma
!  \end{equation}
!  have been used. $\Delta x$ and $\Delta y$ are "small enough", but otherwise
!  arbitrary length scales. The buoyancy gradients computed with this method
!  are then vertically integrated according to \eq{InternalPressurex} and
!  \eq{InternalPressurey}.
!
! The horizontal salinity and temperature gradients have to supplied by the
! user, either as constant values or as profiles given in a file (see
! {\tt obs.nml}).
!
! !USES:
   use meanflow,      only: T,S
   use meanflow,      only: gravity,rho_0,h
   use observations,  only: dsdx,dsdy,dtdx,dtdy
   use observations,  only: idpdx,idpdy,int_press_method
   use eqstate,       only: eqstate1
   IMPLICIT NONE
!
! !INPUT PARAMETERS:

!  number of vertical layers
   integer, intent(in)                 :: nlev
!
! !REVISION HISTORY:
!  Original author(s): Hans Burchard & Karsten Bolding
!
!  $Log: intpressure.F90,v $
!  Revision 1.8  2007-01-06 11:49:15  kbk
!  namelist file extension changed .inp --> .nml
!
!  Revision 1.7  2005/08/11 12:32:50  lars
!  corrected error in Latex referencing
!
!  Revision 1.6  2005/06/27 13:44:07  kbk
!  modified + removed traling blanks
!
!  Revision 1.5  2004/08/18 11:43:51  lars
!  updated documentation
!
!  Revision 1.4  2003/03/28 09:20:35  kbk
!  added new copyright to files
!
!  Revision 1.3  2003/03/28 08:56:56  kbk
!  removed tabs
!
!  Revision 1.2  2003/03/10 08:50:06  gotm
!  Improved documentation and cleaned up code
!
!  Revision 1.1.1.1  2001/02/12 15:55:57  gotm
!  initial import into CVS
!
!EOP
!
! !LOCAL VARIABLES:
   integer                             :: i
   REALTYPE                            :: z,dx,dy
   REALTYPE                            :: dSS,dTT,Bl,Br,int
   REALTYPE                            :: dxB(0:nlev),dyB(0:nlev)
!
!-----------------------------------------------------------------------
!BOC

   if (int_press_method .ne. 0) then

!     initialize local depth
!     and pressure gradient
      z     = _ZERO_
      idpdx = _ZERO_
      idpdy = _ZERO_

!     the spacing for the finite differences
      dx    =  10.
      dy    =  10.

      do i=nlev,1,-1
         z=z+0.5*h(i)

!        buoyancy gradient in x direction
         dSS=dx*dsdx(i)
         dTT=dx*dtdx(i)
         Bl=eqstate1(S(i),T(i),z/10.,gravity,rho_0)
         Br=eqstate1(S(i)+dSS,T(i)+dTT,z/10.,gravity,rho_0)
         dxB(i)=(Br-Bl)/dx

!        buoyancy gradient in y direction
         dSS=dy*dsdy(i)
         dTT=dy*dtdy(i)
         Bl=eqstate1(S(i),T(i),z/10.,gravity,rho_0)
         Br=eqstate1(S(i)+dSS,T(i)+dTT,z/10.,gravity,rho_0)
         dyB(i)=(Br-Bl)/dy

         z=z+0.5*h(i)
      end do

!     internal pressure gradient in x direction
      int=0.5*h(nlev)*dxB(nlev)
      idpdx(nlev)=int
      do i=nlev-1,1,-1
         int=int+0.5*h(i+1)*dxB(i+1)+0.5*h(i)*dxB(i)
         idpdx(i)=int
      end do

!     internal pressure gradient in y direction
      int=0.5*h(nlev)*dyB(nlev)
      idpdy(nlev)=int
      do i=nlev-1,1,-1
         int=int+0.5*h(i+1)*dyB(i+1)+0.5*h(i)*dyB(i)
         idpdy(i)=int
      end do

   endif

   return
   end subroutine intpressure

!EOC

!-----------------------------------------------------------------------
! Copyright by the GOTM-team under the GNU Public License - www.gnu.org
!-----------------------------------------------------------------------