ss_nn.F90 8.51 KB
Newer Older
kbk's avatar
kbk committed
1
!$Id: ss_nn.F90,v 1.7 2006-03-01 15:54:08 kbk Exp $
gotm's avatar
gotm committed
2 3 4 5
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
6
! !ROUTINE: ss_nn - calculates shear and buoyancy frequency
hb's avatar
hb committed
7
! \label{sec-ss-nn}
gotm's avatar
gotm committed
8 9 10 11 12
!
! !INTERFACE:
   subroutine ss_nn()
!
! !DESCRIPTION:
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
!
! Here, the shear frequency squared, 
! $M^2=\left(\partial_z u\right)^2+\left(\partial_z v\right)^2$, 
! and the buoyancy frequency squared, $N^2=\partial_z b$, with buoyancy $b$ from
! (\ref{bdef}) are calculated. 
! For both calculations, two alternative methods are coded.
! The two straight-forward methods which are explained first, do
! both have the disadvantage of generating numerical instabilities.
! The straight-forward way for calculating $M^2$ is as follows:
!
! \begin{equation}\label{ShearSquaredOld}
! \begin{array}{l}
! \displaystyle
! (M^2)_{i,j,k}\approx
! \frac12
! \Bigg(\left(\frac{u_{i,j,k+1}-u_{i,j,k}}
! {\frac12(h^u_{i,j,k+1}+h^u_{i,j,k})}\right)^2
! +
! \left(\frac{u_{i-1,j,k+1}-u_{i-1,j,k}}
! {\frac12(h^u_{i-1,j,k+1}+h^u_{i-1,j,k})}\right)^2
! \\ \\ \qquad\qquad\quad
! \displaystyle
! +
! \left(\frac{v_{i,j,k+1}-v_{i,j,k}}
! {\frac12(h^v_{i,j,k+1}+h^v_{i,j,k})}\right)^2
! +
! \left(\frac{v_{i,j-1,k+1}-v_{i,j-1,k}}
! {\frac12(h^v_{i,j-1,k+1}+h^v_{i,j-1,k})}\right)^2
! \Bigg)
! \end{array}
! \end{equation}
!
! \cite{BURCHARD01c} developed a new scheme, which guarantees that
! the mean kinetic energy which is dissipated from the mean flow
! equals the shear production of turbulent kinetic energy. Therefore,
! this scheme should be numerically more stable than (\ref{ShearSquaredOld}):
!
! \begin{equation}\label{ShearSquaredNew}
! \begin{array}{l}
! \displaystyle
! (M^2)_{i,j,k}\approx
! \frac12
! \Bigg(\frac{\frac12(\nu_{i,j,k}+\nu_{i+1,j,k})
! (u_{i,j,k+1}-u_{i,j,k})^2}{\frac12(h^u_{i,j,k+1}+h^u_{i,j,k})}
! \\ \\ \qquad\qquad\quad
! \displaystyle
! +
! \frac{\frac12(\nu_{i-1,j,k}+\nu_{i,j,k})
! (u_{i-1,j,k+1}-u_{i-1,j,k})^2}{\frac12(h^u_{i-1,j,k+1}+h^u_{i-1,j,k})}
! \\ \\ \qquad\qquad\quad
! \displaystyle
! +
! \frac{\frac12(\nu_{i,j,k}+\nu_{i,j+1,k})
! (v_{i,j,k+1}-v_{i,j,k})^2}{\frac12(h^v_{i,j,k+1}+h^v_{i,j,k})}
! \\ \\ \qquad\qquad\quad
! \displaystyle
! +
! \frac{\frac12(\nu_{i,j-1,k}+\nu_{i,j,k})
! (v_{i,j-1,k+1}-v_{i,j-1,k})^2}{\frac12(h^v_{i,j-1,k+1}+h^v_{i,j-1,k})}
! \Bigg)
! \\ \\ \qquad\qquad\quad
! \displaystyle
! \cdot
! \left(\frac12\left(h^c_{i,j,k}+h^c_{i,j,k+1}\right)\nu_{i,j,k}\right)^{-1}
! \end{array}
! \end{equation}
! 
! The straight-forward discretisation of
! $N^2$ is given by
! 
! \begin{equation}\label{Nstraight}
! \begin{array}{l}
! \displaystyle
! (N^2)_{i,j,k}\approx
! \frac{b_{i,j,k+1}-b_{i,j,k}}{\frac12(h^t_{i,j,k+1}+h^t_{i,j,k})}.
! \end{array}
! \end{equation}
! 
! In some cases, together with the straight-forward discretisation
! of the shear squared, (\ref{ShearSquaredOld}), this
! did not produce stable numerical results. The reason for this might be that
! the velocities involved in the calculation for the shear squared do depend
! on the buoyancies in the two
! neighbouring T-points such that the straight-forward
! method (\ref{Nstraight}) leads to an inconsistency.
! However, other experiments with the energy-conserving discretisation
! of the shear stress squared, (\ref{ShearSquaredNew})
! and the straight-forward discretisation of
! $N^2$, (\ref{Nstraight}),  produced numerically stable results.
!
! Most stable results have been obtained with a weighted average 
! for the $N^2$ calculation:
!
! \begin{equation}\label{Naveraged}
! \begin{array}{l}
! \displaystyle
! (N^2)_{i,j,k}\approx
! \frac16 \Bigg(
! 2\frac{b_{i,j,k+1}-b_{i,j,k}}{\frac12(h^t_{i,j,k+1}+h^t_{i,j,k})}
! \\ \\ \qquad\qquad\qquad
! \displaystyle
! +
! \frac{b_{i+1,j,k+1}-b_{i+1,j,k}}{\frac12(h^t_{i+1,j,k+1}+h^t_{i+1,j,k})}
! +
! \frac{b_{i-1,j,k+1}-b_{i-1,j,k}}{\frac12(h^t_{i-1,j,k+1}+h^t_{i-1,j,k})}
! \\ \\ \qquad\qquad\qquad
! \displaystyle
! +
! \frac{b_{i,j+1,k+1}-b_{i,j+1,k}}{\frac12(h^t_{i,j+1,k+1}+h^t_{i,j+1,k})}
! +
! \frac{b_{i,j-1,k+1}-b_{i,j-1,k}}{\frac12(h^t_{i,j-1,k+1}+h^t_{i,j-1,k})}
! \Bigg).
! \end{array}
! \end{equation}
!
! These stability issues need to be further investigated in the future. 
gotm's avatar
gotm committed
129 130
!
! !USES:
131
   use domain, only: iimin,iimax,jjmin,jjmax,kmax,au,av,az
132
   use variables_3d, only: kmin,kumin,hn,uu,hun,kvmin,vv,hvn,SS,num
133
#ifndef NO_BAROCLINIC
134
   use variables_3d, only: NN,rho
135
#endif
gotm's avatar
gotm committed
136 137 138 139 140 141 142 143
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
kbk's avatar
kbk committed
144 145 146
! !REVISION HISTORY:
!  Original author(s): Hans Burchard & Karsten Bolding
!
gotm's avatar
gotm committed
147
! !LOCAL VARIABLES:
kbk's avatar
kbk committed
148 149
   integer                   :: i,j,k,nb
   REALTYPE                  :: dz,NNc,NNe,NNw,NNn,NNs
gotm's avatar
gotm committed
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
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'ss_nn() # ',Ncall
#endif

#undef NEW_SS
#define NEW_SS

!  Prandtl frequency
   do j=jjmin,jjmax
      do i=iimin,iimax
         if (az(i,j) .eq. 1 ) then
            do k=1,kmax-1
! This is an older version which we should keep here.
#ifndef NEW_SS
              SS(i,j,k)=0.5* (                                               &
                   ( (uu(i,j,k+1)/hun(i,j,k+1)-uu(i,j,k)/hun(i,j,k))         &
                   /(0.5*(hun(i,j,k+1)+hun(i,j,k))) )**2                     &
                +  ( (uu(i-1,j,k+1)/hun(i-1,j,k+1)-uu(i-1,j,k)/hun(i-1,j,k)) &
                   /(0.5*(hun(i-1,j,k+1)+hun(i-1,j,k))) )**2                 &
                +  ( (vv(i,j,k+1)/hvn(i,j,k+1)-vv(i,j,k)/hvn(i,j,k))         &
                   /(0.5*(hvn(i,j,k+1)+hvn(i,j,k))) )**2                     &
                +  ( (vv(i,j-1,k+1)/hvn(i,j-1,k+1)-vv(i,j-1,k)/hvn(i,j-1,k)) &
                   /(0.5*(hvn(i,j-1,k+1)+hvn(i,j-1,k))) )**2                 &
                            )
#else
! This version should better conserve energy.
              SS(i,j,k)=0.5* (                                                &
                   (uu(i,j,k+1)/hun(i,j,k+1)-uu(i,j,k)/hun(i,j,k))**2         &
                   /(0.5*(hun(i,j,k+1)+hun(i,j,k)))                           &
kbk's avatar
kbk committed
184
                    *0.5*(num(i,j,k)+num(i+1,j,k))                            &
gotm's avatar
gotm committed
185 186 187 188 189
               +  (uu(i-1,j,k+1)/hun(i-1,j,k+1)-uu(i-1,j,k)/hun(i-1,j,k))**2 &
                  /(0.5*(hun(i-1,j,k+1)+hun(i-1,j,k)))                       &
                   *0.5*(num(i-1,j,k)+num(i,j,k))                            &
                +  (vv(i,j,k+1)/hvn(i,j,k+1)-vv(i,j,k)/hvn(i,j,k))**2         &
                   /(0.5*(hvn(i,j,k+1)+hvn(i,j,k)))                           &
kbk's avatar
kbk committed
190
                    *0.5*(num(i,j,k)+num(i,j+1,k))                            &
gotm's avatar
gotm committed
191 192
                +  (vv(i,j-1,k+1)/hvn(i,j-1,k+1)-vv(i,j-1,k)/hvn(i,j-1,k))**2 &
                   /(0.5*(hvn(i,j-1,k+1)+hvn(i,j-1,k)))                       &
kbk's avatar
kbk committed
193
                    *0.5*(num(i,j-1,k)+num(i,j,k))                            &
gotm's avatar
gotm committed
194 195 196 197 198 199 200
                            )/(0.5*(hn(i,j,k)+hn(i,j,k+1)))/num(i,j,k)
#endif
            end do
         end if
      end do
   end do

201
#ifndef NO_BAROCLINIC
gotm's avatar
gotm committed
202 203 204 205 206 207
#define NEW_NN
#undef NEW_NN
   do j=jjmin,jjmax
      do i=iimin,iimax
         if (az(i,j) .eq. 1 ) then
            do k=kmax-1,1,-1
kbk's avatar
kbk committed
208
               dz=0.5*(hn(i,j,k+1)+hn(i,j,k))
gotm's avatar
gotm committed
209 210
               NNc =(rho(i,j,k+1)-rho(i,j,k))/dz
#ifndef NEW_NN
kbk's avatar
kbk committed
211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235
               if (az(i+1,j) .eq. 1) then
                  dz=0.5*(hn(i+1,j,k+1)+hn(i+1,j,k))
                  NNe=(rho(i+1,j,k+1)-rho(i+1,j,k))/dz
               else
                  NNe=NNc
               end if 
               if (az(i-1,j) .eq. 1) then
                  dz=0.5*(hn(i-1,j,k+1)+hn(i-1,j,k))
                  NNw=(rho(i-1,j,k+1)-rho(i-1,j,k))/dz
               else
                  NNw=NNc
               end if
               if (az(i,j+1) .eq. 1) then
                  dz=0.5*(hn(i,j+1,k+1)+hn(i,j+1,k))
                  NNn=(rho(i,j+1,k+1)-rho(i,j+1,k))/dz
               else
                  NNn=NNc
               end if
               if (az(i,j-1) .eq. 1) then
                  dz=0.5*(hn(i,j-1,k+1)+hn(i,j-1,k))
                  NNs=(rho(i,j-1,k+1)-rho(i,j-1,k))/dz
               else
                  NNs=NNc
               end if
               NN(i,j,k)=0.3333333*NNc+0.1666666*(NNe+NNw+NNn+NNs)
gotm's avatar
gotm committed
236 237 238 239 240 241 242
#else
               NN(i,j,k)=NNc
#endif
            end do
         end if
      end do
   end do
243
#endif
gotm's avatar
gotm committed
244 245 246 247 248 249 250 251 252 253 254 255

#ifdef DEBUG
   write(debug,*) 'Leaving ss_nn()'
   write(debug,*)
#endif
   return
   end subroutine ss_nn
!EOC

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