ip_song_wright.F90 4.54 KB
Newer Older
kbk's avatar
kbk committed
1
!$Id: ip_song_wright.F90,v 1.4 2006-03-01 15:54:08 kbk Exp $
2 3 4 5
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
6
! !ROUTINE: ip_song_wright
7 8 9 10 11 12
!
! !INTERFACE:
   subroutine ip_song_wright()
!
! !DESCRIPTION:
!
hb's avatar
hb committed
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
! Here, the pressure gradient is calculating according to an energy-conserving
! method suggested by \cite{SONG98}, which for the pressure gradient in
! $x$-direction looks as:
! 
! \begin{equation}\label{drhodxdiscrSONG}
! \begin{array}{l}
! \displaystyle
! \frac12(h_{i,j,k}+h_{i,j,k+1})\left(m\,\partial_{\cal X}^*b\right)_{i,j,k} \\ \\
! \displaystyle
! \approx
! \frac{
! \frac14 (b_{i+1,j,k+1}+b_{i+1,j,k})(h^c_{i+1,j,k+1}+h^c_{i+1,j,k})-
! \frac14 (b_{i,j,k+1}+b_{i,j,k})(h^c_{i,j,k+1}+h^c_{i,j,k})}
! {\Delta x^u_{i,j}}\\ \\
! \displaystyle
! \qquad-
! \Bigg[\frac12 (b_{i+1,j,k+1}+b_{i,j,k+1})
! \frac{z^c_{i+1,j,k+1}-z^c_{i,j,k+1}}{\Delta x^u_{i,j}}\\ \\
! \displaystyle
! \qquad\qquad -
! \frac12 (b_{i+1,j,k}+b_{i,j,k})
! \frac{z^c_{i+1,j,k}-z^c_{i,j,k}}{\Delta x^u_{i,j}}\Bigg],
! \end{array}
! \end{equation}
! 
! where $z^c_{i,j,k}$ is the $z$-coordinate of the centre of
! the grid box with the index $(i,j,k)$.
!
! The discretisation of $(\partial_y^* b)_k$ for the $v$-equation is
! done accordingly.
!
!
45 46 47 48
! !USES:
   use internal_pressure
   IMPLICIT NONE
!
kbk's avatar
kbk committed
49 50 51
! !REVISION HISTORY:
!  Original author(s): Hans Burchard & Karsten Bolding
!
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
! !LOCAL VARIABLES:
   integer                   :: i,j,k
   REALTYPE                  :: dxm1,dym1
   REALTYPE                  :: grdl,grdu,rhol,rhou,prgr,dxz,dyz
!EOP
!-----------------------------------------------------------------------
!BOC
#if ! ( defined(SPHERICAL) || defined(CURVILINEAR) )
   dxm1 = _ONE_/DXU
   dym1 = _ONE_/DYV
#endif

!  First, the interface heights are calculated in order to get the
!  interface slopes further down.
   do j=jjmin,jjmax+1
      do i=iimin,iimax+1
         if (az(i,j) .ge. 1) then
            zz(i,j,1)=-H(i,j)+0.5*hn(i,j,1)
            do k=2,kmax
               zz(i,j,k)=zz(i,j,k-1)+0.5*(hn(i,j,k-1)+hn(i,j,k))
            end do
         end if
      end do
   end do

!  Calculation of layer integrated internal pressure gradient as it
!  appears on the right hand side of the u-velocity equation.
   do j=jjmin,jjmax
      do i=iimin,iimax
         if (au(i,j) .ge. 1) then
#if defined(SPHERICAL) || defined(CURVILINEAR)
            dxm1=_ONE_/DXU
#endif
            grdl=0.5*hun(i,j,kmax)*(rho(i+1,j,kmax)-rho(i,j,kmax))*dxm1
            rhol=0.5*(rho(i+1,j,kmax)+rho(i,j,kmax))       &
                     *(zz(i+1,j,kmax)- zz(i,j,kmax))*dxm1
            prgr=grdl
            idpdx(i,j,kmax)=hun(i,j,kmax)*prgr
            do k=kmax-1,1,-1
               grdu=grdl
               grdl=(0.5*(rho(i+1,j,k)+rho(i+1,j,k+1))               &
                        *0.5*(hn(i+1,j,k)+hn(i+1,j,k+1))             &
                    -0.5*(rho(i,j,k)+rho(i,j,k+1))                   &
                        *0.5*(hn(i,j,k)+hn(i,j,k+1)) )*dxm1
               rhou=rhol
               rhol=0.5*(rho(i+1,j,k)+rho(i,j,k))*(zz(i+1,j,k)-zz(i,j,k))*dxm1
               prgr=prgr+grdl-(rhou-rhol)
               idpdx(i,j,k)=hun(i,j,k)*prgr
            end do
         end if
      end do
   end do

!  Calculation of layer integrated internal pressure gradient as it
!  appears on the right hand side of the v-velocity equation.
   do j=jjmin,jjmax
      do i=iimin,iimax
         if (av(i,j) .ge. 1) then
#if defined(SPHERICAL) || defined(CURVILINEAR)
         dym1 = _ONE_/DYV
#endif
            grdl=0.5*hvn(i,j,kmax)*(rho(i,j+1,kmax)-rho(i,j,kmax))*dym1
            rhol=0.5*(rho(i,j+1,kmax)+rho(i,j,kmax))     &
                     *(zz(i,j+1,kmax)- zz(i,j,kmax))*dxm1
            prgr=grdl
            idpdy(i,j,kmax)=hvn(i,j,kmax)*prgr
            do k=kmax-1,1,-1
               grdu=grdl
               grdl=(0.5*(rho(i,j+1,k)+rho(i,j+1,k+1))               &
                        *0.5*(hn(i,j+1,k)+hn(i,j+1,k+1))             &
                    -0.5*(rho(i,j,k)+rho(i,j,k+1))                   &
                        *0.5*(hn(i,j,k)+hn(i,j,k+1)) )*dym1
               rhou=rhol
               rhol=0.5*(rho(i,j+1,k)+rho(i,j,k))*(zz(i,j+1,k)-zz(i,j,k))*dym1
               prgr=prgr+grdl-(rhou-rhol)
               idpdy(i,j,k)=hvn(i,j,k)*prgr
            end do
         end if
      end do
   end do

   return
   end subroutine ip_song_wright
!EOC

!-----------------------------------------------------------------------
! Copyright (C) 2004 - Hans Burchard and Karsten Bolding               !
!-----------------------------------------------------------------------