extpressure.F90 4.58 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
!$Id: extpressure.F90,v 1.9 2007-01-06 11:57:08 kbk Exp $
#include"cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: The external pressure-gradient \label{sec:extpressure}
!
! !INTERFACE:
   subroutine extpressure(method,nlev)
!
! !DESCRIPTION:
!
!  This subroutine calculates the external pressure-gradient. Two methods
!  are implemented here, relating either to the velocity vector at a
!  given height above bed prescribed or to the vector for the vertical mean
!  velocity. In the first case, {\tt dpdx} and {\tt dpdy} are $x$-
!  and $y$-components of the prescribed velocity vector at the
!  height {\tt h\_press} above the bed. The velocity profile will in
!  this routive be shifted by a vertically constant vector such that the
!  resulting profile has an (interpolated) velocity at {\tt h\_press}
!  which is identical to the prescribed value. In the second case,
!  {\tt dpdx} and {\tt dpdy} are $x$- and $y$-components of the
!  prescribed vertical mean velocity vector, and {\tt h\_press} is
!  not used. Here the velocity profile is shifted in such a way that
!  the resulting mean velocty vector is identical to {\tt dpdx} and {\tt dpdy}.
!
!  For both cases, this is a recalculation of the external pressure gradient,
!  since at all points the same acceleration has been applied in this
!  operator split method.
!
!  If the external pressure-gradient is prescribed by the
!  surface slope, then it is directly inserted in \eq{uEq} and \eq{vEq}.
!
!  For details of this method, see \cite{Burchard99}.
!
! !USES:
   use meanflow,     only: u,v,h
   use observations, only: dpdx,dpdy,h_press
!
   IMPLICIT NONE
!
! !INPUT PARAMETERS:

!  method to compute external
!  pressure gradient
   integer, intent(in)                 :: method

!  number of vertical layers
   integer, intent(in)                 :: nlev
!
! !REVISION HISTORY:
!  Original author(s): Hans Burchard & Karsten Bolding
!
!  $Log: extpressure.F90,v $
!  Revision 1.9  2007-01-06 11:57:08  kbk
!  PressMethod --> ext_press_mode
!
!  Revision 1.8  2007-01-06 11:49:15  kbk
!  namelist file extension changed .inp --> .nml
!
!  Revision 1.7  2006-12-07 16:50:28  hb
!  Program abortion introduced for PressHeight<=0 when ext_press_mode=1
!
!  Revision 1.6  2005-06-27 13:44:07  kbk
!  modified + removed traling blanks
!
!  Revision 1.5  2004/08/18 11:41:02  lars
!  corrected typo in docu
!
!  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(0:nlev)
   REALTYPE                            :: rat,uint,vint,hint
!
!-----------------------------------------------------------------------
!BOC

   select case (method)
      case (1)
!        current measurement at h_press above bed
         if (h_press.le._ZERO_) then
            LEVEL2 ''
            LEVEL2 '***************************************************'
            LEVEL2 'PressHeight=',h_press, 'but it must be positive.'
            LEVEL2 'Please correct this in obs.nml and rerun.'
            LEVEL2 'Program aborted.'
            LEVEL2 '***************************************************' 
            stop 'extpressure'
         end if
         z(1)=0.5*h(1)
         i   =0
222      i=i+1
         z(i+1)=z(i)+0.5*(h(i)+h(i+1))
         if ((z(i+1).lt.h_press).and.(i.lt.nlev)) goto 222
         rat=(h_press-z(i))/(z(i+1)-z(i))
         uint=rat*u(i+1)+(1-rat)*u(i)
         vint=rat*v(i+1)+(1-rat)*v(i)
         do i=1,nlev
            u(i)=u(i)+dpdx-uint
            v(i)=v(i)+dpdy-vint
         end do
      case (2)
!     vertical mean of current prescribed
         uint=0.
         vint=0.
         hint=0.
         do i=1,nlev
            hint=hint+h(i)
            uint=uint+h(i)*u(i)
            vint=vint+h(i)*v(i)
         end do
         uint=uint/hint
         vint=vint/hint
         do i=1,nlev
            u(i)=u(i)+dpdx-uint
            v(i)=v(i)+dpdy-vint
         end do
      case default
!     do nothing if method=0, because then
!     pressure gradient is applied directly
!     in uequation() and vequation()
   end select

   return
   end subroutine extpressure
!EOC

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