buoyancy.F90 5.31 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
!$Id: buoyancy.F90,v 1.8 2006-11-06 13:36:45 hb Exp $
#include"cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: The buoyancy equation
!
! !INTERFACE:
   subroutine buoyancy(nlev,dt,cnpar,nub,gamb)
!
! !DESCRIPTION:
!  This subroutine solves a transport equation for the mean
!  potential buoyancy,
!  \begin{equation}
!    \label{DefBuoyancy}
!    B=-g\frac{\mean{\rho}-\rho_0}{\rho_0}
!    \comma
!  \end{equation}
!  where $g$ is the accelaration of gravity, and $\mean{\rho}$ and $\rho_0$
!  are the mean potential density and the reference density, respectively.
!  A simplified transport equation for $B$ can be written as
!  \begin{equation}
!   \label{bEq}
!    \dot{B}
!    = {\cal D}_B
!    \comma
!  \end{equation}
!  where $\dot{B}$ denotes the material derivative of $B$, and
!  ${\cal D}_b$ is the sum of the turbulent and viscous transport
!  terms modelled according to
!  \begin{equation}
!   \label{Db}
!    {\cal D}_B
!    = \frstder{z}
!     \left( (\nu^B_t+\nu^B) \partder{B}{z} - \tilde{\Gamma}_B \right)
!    \point
!  \end{equation}
!  In this equation, $\nu^B_t$ and $\nu^B$ are the turbulent and molecular
!  diffusivities of buoyancy, respectively, and $\tilde{\Gamma}_B$
!  denotes the non-local flux of buoyancy, see
!  \sect{sec:turbulenceIntro}. In the current version of GOTM,
!  we set $\nu^B_t = \nu^\Theta_t$ for simplicity. Source and sink
!  terms are completely disregarded, and thus \eq{bEq} mainly serves
!  as a convenient tool for some idealized test cases in
!  GOTM.
!
!  Diffusion is treated implicitly in space (see equations (\ref{sigmafirst})-
!  (\ref{sigmalast})), and then solved by a
!  simplified Gauss elimination.
!  Vertical advection is included, and it must be non-conservative,
!  which is ensured by setting the local variable {\tt adv\_mode=0},
!  see section \ref{sec:advectionMean} on page \pageref{sec:advectionMean}.
!
! !USES:
   use meanflow,      only: h,w,buoy,T,avh
   use meanflow,      only: w_grid,grid_method
   use observations,  only: b_obs_NN,b_obs_surf,b_obs_sbf
   use observations,  only: w_adv_discr,w_adv_method
   use util,          only: Dirichlet,Neumann
   use util,          only: oneSided,zeroDivergence
!
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
!  number of vertical layers
   integer, intent(in)                 :: nlev

!   time step (s)
   REALTYPE, intent(in)                :: dt

!  numerical "implicitness" parameter
   REALTYPE, intent(in)                :: cnpar

!  diffusivity of buoyancy (m^2/s)
   REALTYPE, intent(in)                :: nub(0:nlev)

!  non-local buoyancy flux (m^2/s^3)
   REALTYPE, intent(in)                :: gamb(0:nlev)
!
!
! !REVISION HISTORY:
!  Original author(s): Hans Burchard & Karsten Bolding
!
!  $Log: buoyancy.F90,v $
!  Revision 1.8  2006-11-06 13:36:45  hb
!  Option for conservative vertical advection added to adv_center
!
!  Revision 1.7  2005-11-17 09:58:20  hb
!  explicit argument for positive definite variables in diff_center()
!
!  Revision 1.6  2005/06/27 13:44:07  kbk
!  modified + removed traling blanks
!
!  Revision 1.5  2003/03/28 09:20:35  kbk
!  added new copyright to files
!
!  Revision 1.4  2003/03/28 08:56:56  kbk
!  removed tabs
!
!  Revision 1.3  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                   :: adv_mode=0
   integer                   :: posconc=0
   integer                   :: i
   integer                   :: DiffBcup,DiffBcdw
   integer                   :: AdvBcup,AdvBcdw
   REALTYPE                  :: DiffBup,DiffBdw
   REALTYPE                  :: AdvBup,AdvBdw
   REALTYPE                  :: Lsour(0:nlev)
   REALTYPE                  :: Qsour(0:nlev)
   REALTYPE                  :: BRelaxTau(0:nlev)
   REALTYPE                  :: zz

   logical, save             :: first=.true.
!
!-----------------------------------------------------------------------
!BOC

!  no diffusive flux at bottom and surface
   DiffBcup      = Neumann
   DiffBcdw      = Neumann
   DiffBup       = _ZERO_
   DiffBdw       = _ZERO_

   AdvBcup       = zeroDivergence
   AdvBcdw       = oneSided
   AdvBup        = _ZERO_
   AdvBdw        = _ZERO_


!  construct initial linear profile from information in namelist
   if (first) then

      zz=_ZERO_
      do i=nlev,1,-1
         zz=zz+0.5*h(i)
         buoy(i)  = b_obs_surf - zz*b_obs_NN
         zz=zz+0.5*h(i)
      end do

      first=.false.
   end if

!  compose source term
   do i=1,nlev
      Lsour(i) = _ZERO_
      Qsour(i) = - (gamB(i)-gamB(i-1))/h(i)
   end do

!  no relaxation to observed values
   BRelaxTau = 1.e15

!  do advection step
   if (w_adv_method .ne. 0) then
      call adv_center(nlev,dt,h,h,w,AdvBcup,AdvBcdw,                    &
                      AdvBup,AdvBdw,w_adv_discr,adv_mode,buoy)
   end if

!  do diffusion step
   call diff_center(nlev,dt,cnpar,posconc,h,DiffBcup,DiffBcdw,          &
                    DiffBup,DiffBdw,avh,Lsour,Qsour,                    &
                    BRelaxTau,buoy,buoy)

!   T = buoy

   return
   end subroutine buoyancy
!EOC

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