uv_diffusion.F90 9.5 KB
Newer Older
kbk's avatar
kbk committed
1
!$Id: uv_diffusion.F90,v 1.8 2006-03-01 15:54:07 kbk Exp $
gotm's avatar
gotm committed
2 3 4 5
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
6
! !ROUTINE: uv_diffusion - 2D diffusion of momentum \label{sec-uv-diffusion}
gotm's avatar
gotm committed
7 8
!
! !INTERFACE:
9
   subroutine uv_diffusion(Am,An)
gotm's avatar
gotm committed
10 11 12
!
! !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 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
! Here, the diffusion terms for the vertically integrated transports are
! calculated by means of central differences, following the finite volume
! approach. They are added to the
! advection terms into the terms {\tt UEx} and {\tt VEx} for the
! $U$- and the $V$-equation, respectively. The physical diffusion with the
! given eddy viscosity coefficient $A_h^M$ is based on velocity gradients,
! whereas an additional numerical damping of the barotropic mode is based
! on gradients of the transports with the damping coefficient $A_h^N$, 
! see the example given as equations (\ref{smooth_example_1}) and 
! (\ref{smooth_example_2}).
!
! First diffusion term in (\ref{UMOM}):
! \begin{equation}
! \left(mn\,\partial_{\cal X}\left(2A_h^MD\partial_{\cal X}\left(\frac{U}{D}\right)
! + A_h^N\partial_{\cal X}U\right)\right)_{i,j}\approx 
! \frac{
! {\cal F}^{Dxx}_{i+1,j}-{\cal F}^{Dxx}_{i,j}
! }{\Delta x^u_{i,j}\Delta y^u_{i,j}}
! \end{equation}
!
! with diffusive fluxes
!
! \begin{equation}
! {\cal F}^{Dxx}_{i,j}=\left(2 A_h^MD_{i,j}\left(\frac{U_{i,j}}{D^u_{i,j}}
! -\frac{U_{i-1,j}}{D^u_{i-1,j}}\right)+A_h^N\left(U_{i,j}
! -U_{i-1,j}\right)\right)
! \frac{\Delta y^c_{i,j}}{\Delta x^c_{i,j}}.
! \end{equation}
!
! Second diffusion term in (\ref{UMOM}):
! \begin{equation}
! \left(mn\,\partial_{\cal Y}\left(A_h^MD\left(\partial_{\cal Y}\left(\frac{U}{D}\right)+\partial_{\cal X}\left(\frac{V}{D}\right)\right)
! + A_h^N\partial_{\cal Y}U\right)\right)_{i,j}\approx 
! \frac{
! {\cal F}^{Dxy}_{i,j}-{\cal F}^{Dxy}_{i,j-1}
! }{\Delta x^x_{i,j}\Delta y^x_{i,j}}
! \end{equation}
!
! with diffusive fluxes
!
! \begin{equation}
! \begin{array}{rcl}
! \displaystyle
! {\cal F}^{Dxy}_{i,j}&=&
! \displaystyle
! A_h^M\frac12\left(D^u_{i,j}+D^u_{i,j+1}\right)
! \Delta x^x_{i,j}\left(\left(\frac{U_{i,j+1}}{D^u_{i,j+1}}
! -\frac{U_{i,j}}{D^u_{i,j}}\right)\frac{1}{\Delta y^x_{i,j}}
! +\left(\frac{V_{i+1,j}}{D^v_{i+1,j}}
! -\frac{V_{i,j}}{D^v_{i,j}}\right)\frac{1}{\Delta x^x_{i,j}}\right) \\ \\
! &&
! \displaystyle
! +A_h^N\left(U_{i,j+1} -U_{i,j}\right)\frac{\Delta x^x_{i,j}}
! {\Delta y^x_{i,j}}.
! \end{array}
! \end{equation}
!
! First diffusion term in (\ref{VMOM}):
! \begin{equation}
! \left(mn\,\partial_{\cal X}\left(A_h^MD\left(\partial_{\cal Y}\left(\frac{U}{D}\right)+\partial_{\cal X}\left(\frac{V}{D}\right)\right)
! + A_h^N\partial_{\cal X}V\right)\right)_{i,j}\approx 
! \frac{
! {\cal F}^{Dyx}_{i,j}-{\cal F}^{Dyx}_{i-1,j}
! }{\Delta x^x_{i,j}\Delta y^x_{i,j}}
! \end{equation}
!
! with diffusive fluxes
!
! \begin{equation}
! \begin{array}{rcl}
! \displaystyle
! {\cal F}^{Dyx}_{i,j}&=&
! \displaystyle
! A_h^M\frac12\left(D^v_{i,j}+D^v_{i+1,j}\right)
! \Delta y^x_{i,j}\left(\left(\frac{U_{i,j+1}}{D^u_{i,j+1}}
! -\frac{U_{i,j}}{D^u_{i,j}}\right)\frac{1}{\Delta y^x_{i,j}}
! +\left(\frac{V_{i+1,j}}{D^v_{i+1,j}}
! -\frac{V_{i,j}}{D^v_{i,j}}\right)\frac{1}{\Delta x^x_{i,j}}\right) \\ \\
! &&
! \displaystyle
! +A_h^N\left(V_{i+1,j} -V_{i,j}\right)\frac{\Delta y^x_{i,j}}
! {\Delta x^x_{i,j}}.
! \end{array}
! \end{equation}
!
! Second diffusion term in (\ref{VMOM}):
! \begin{equation}
! \left(mn\,\partial_{\cal Y}\left(2A_h^MD\partial_{\cal Y}\left(\frac{V}{D}\right)
! + A_h^N\partial_{\cal Y}V\right)\right)_{i,j}\approx 
! \frac{
! {\cal F}^{Dyy}_{i,j+1}-{\cal F}^{Dyy}_{i,j}
! }{\Delta x^v_{i,j}\Delta y^v_{i,j}}
! \end{equation}
!
! with diffusive fluxes
!
! \begin{equation}
! {\cal F}^{Dyy}_{i,j}=\left(2 A_h^MD_{i,j}\left(\frac{V_{i,j}}{D^v_{i,j}}
! -\frac{V_{i,j-1}}{D^v_{i,j-1}}\right)+A_h^N\left(V_{i,j}
! -V_{i,j-1}\right)\right)
! \frac{\Delta x^c_{i,j}}{\Delta y^c_{i,j}}.
! \end{equation}
!
! The role of the additional diffusion of $U$ and $V$ with the
! diffusion coefficient $A_h^N$ is best demonstrated by means of a
! simplified set of vertically integrated equations:
!
! \begin{equation}\label{smooth_example_1}
! \begin{array}{l}
! \displaystyle
! \partial_t \eta = - \partial_x U - \partial_y V \\ \\
! \displaystyle
! \partial_t U = -gD\partial_x\eta 
! + A_h^N \left(\partial_{xx} U + \partial_{yy} U\right) \\ \\
! \displaystyle
! \partial_t V = -gD\partial_y\eta 
! + A_h^N \left(\partial_{xx} V + \partial_{yy} V\right), \\ \\
! \end{array}
! \end{equation}
!
! which can be transformed into an equation for $\partial_t\eta$ by
! derivation of the $\eta$-equation with respect to $t$, 
! of the $U$-equation with respect to $x$
! and the $V$-equation with respect to $y$ and subsequent elimination of 
! $U$ and $V$:
!
! \begin{equation}\label{smooth_example_2}
! \partial_t \left(\partial_t\eta\right) = gD \left(\partial_{xx}\eta
! + \partial_{yy}\eta\right)
! +  A_h^N \left(\partial_{xx}\left(\partial_t\eta\right)
! +\partial_{yy}\left(\partial_t\eta\right)\right),
! \end{equation}
!
! which can be interpreted as a wave equation with a damping on
! $\partial_t \eta$. This introduces an explicit damping of
! free surface elevation oscillations in a momentum-conservative
! manner. Hydrodynamic models with implicit treatment of the barotropic mode
! do not need to apply this method due to the implicit damping of those models,
! see e.g.\ \cite{BACKHAUS85}. The implementation of this explicit damping 
! described here has been suggested by Jean-Marie Beckers, Li\'ege
! (Belgium).
! 
! When working with the option {\tt SLICE\_MODEL}, the calculation of
! all gradients in $y$-direction is suppressed.
!
gotm's avatar
gotm committed
158
! !USES:
159
   use domain, only: imin,imax,jmin,jmax,az,au,av,ax
gotm's avatar
gotm committed
160 161 162 163 164
#if defined(SPHERICAL) || defined(CURVILINEAR)
   use domain, only: dyc,arud1,dxx,dyx,arvd1,dxc
#else
   use domain, only: dx,dy,ard1
#endif
kbk's avatar
kbk committed
165
   use variables_2d, only: D,U,DU,UEx,V,DV,VEx,PP
gotm's avatar
gotm committed
166 167 168
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
169
  REALTYPE, intent(in) :: Am,An
gotm's avatar
gotm committed
170 171 172 173 174
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
kbk's avatar
kbk committed
175 176 177
! !REVISION HISTORY:
!  Original author(s): Hans Burchard
!
gotm's avatar
gotm committed
178
! !LOCAL VARIABLES:
179
   integer                   :: i,j
gotm's avatar
gotm committed
180 181 182 183 184 185 186 187 188
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'uv_diffusion() # ',Ncall
#endif

189
! Central for dx(2*Am*dx(U/DU))
gotm's avatar
gotm committed
190 191
   do j=jmin,jmax
      do i=imin,imax+1          ! PP defined on T-points
192
         PP(i,j)=_ZERO_
gotm's avatar
gotm committed
193
         if (az(i,j) .ge. 1) then
194 195 196 197 198 199 200
            if(Am .gt. _ZERO_) then
               PP(i,j)=2.*Am*DYC*D(i,j)               &
                       *(U(i,j)/DU(i,j)-U(i-1,j)/DU(i-1,j))/DXC
            end if
            if(An .gt. _ZERO_) then
               PP(i,j)=PP(i,j)+An*DYC*(U(i,j)-U(i-1,j))/DXC
            end if
gotm's avatar
gotm committed
201 202 203 204 205 206 207
         end if
      end do
   end do
   do j=jmin,jmax      ! UEx defined on U-points
      do i=imin,imax
         if (au(i,j) .ge. 1) then
            UEx(i,j)=UEx(i,j)-(PP(i+1,j)-PP(i  ,j))*ARUD1
kbk's avatar
kbk committed
208
         end if
gotm's avatar
gotm committed
209 210 211
      end do
   end do

212
#ifndef SLICE_MODEL
213
! Central for dy(Am*(dy(U/DU)+dx(V/DV)))
gotm's avatar
gotm committed
214
   do j=jmin-1,jmax        ! PP defined on X-points
215
      do i=imin,imax
216
         PP(i,j)=_ZERO_
217
         if (ax(i,j) .ge. 1) then
218 219 220 221 222 223
            if(Am .gt. _ZERO_) then
               PP(i,j)=Am*0.5*(DU(i,j)+DU(i,j+1))*DXX  &
                       *((U(i,j+1)/DU(i,j+1)-U(i,j)/DU(i,j))/DYX &
                        +(V(i+1,j)/DV(i+1,j)-V(i,j)/DV(i,j))/DXX )
            end if
            if(An .gt. _ZERO_) then
224
               PP(i,j)=PP(i,j)+An*(U(i,j+1)-U(i,j))*DXX/DYX
225
            end if
kbk's avatar
kbk committed
226
         end if
gotm's avatar
gotm committed
227 228 229 230 231 232
      end do
   end do
   do j=jmin,jmax        !UEx defined on U-points
      do i=imin,imax
         if (au(i,j) .ge. 1) then
            UEx(i,j)=UEx(i,j)-(PP(i,j  )-PP(i,j-1))*ARUD1
kbk's avatar
kbk committed
233
         end if
gotm's avatar
gotm committed
234 235
      end do
   end do
236
#endif
gotm's avatar
gotm committed
237 238

! Central for dx(Am*(dy(U^2/DU)+dx(V^2/DV)))
239
   do j=jmin,jmax      ! PP defined on X-points
gotm's avatar
gotm committed
240
      do i=imin-1,imax
241
         PP(i,j)=_ZERO_
242
         if (ax(i,j) .ge. 1) then
243 244 245 246 247 248
            if(Am .gt. _ZERO_) then
               PP(i,j)=Am*0.5*(DV(i,j)+DV(i+1,j))*DYX  &
                       *((U(i,j+1)/DU(i,j+1)-U(i,j)/DU(i,j))/DYX &
                        +(V(i+1,j)/DV(i+1,j)-V(i,j)/DV(i,j))/DXX )
            end if
            if(An .gt. _ZERO_) then
249
               PP(i,j)=PP(i,j)+An*(V(i+1,j)-V(i,j))*DYX/DXX
250
            end if
kbk's avatar
kbk committed
251
         end if
gotm's avatar
gotm committed
252 253 254 255 256 257
      end do
   end do
   do j=jmin,jmax          ! VEx defined on V-points
      do i=imin,imax
         if (av(i,j) .ge. 1) then
            VEx(i,j)=VEx(i,j)-(PP(i  ,j)-PP(i-1,j))*ARVD1
kbk's avatar
kbk committed
258
         end if
gotm's avatar
gotm committed
259 260 261
      end do
   end do

262
#ifndef SLICE_MODEL
gotm's avatar
gotm committed
263 264
! Central for dy(2*Am*dy(V^2/DV))
   do j=jmin,jmax+1     ! PP defined on T-points
265
      do i=imin,imax
266
         PP(i,j)=_ZERO_
kbk's avatar
kbk committed
267
         if (az(i,j) .ge. 1) then
268 269 270 271 272 273 274
            if(Am .gt. _ZERO_) then
               PP(i,j)=2.*Am*DXC*D(i,j)               &
                       *(V(i,j)/DV(i,j)-V(i,j-1)/DV(i,j-1))/DYC
            end if
            if(An .gt. _ZERO_) then
               PP(i,j)=PP(i,j)+An*DXC*(V(i,j)-V(i,j-1))/DYC
            end if
kbk's avatar
kbk committed
275
         end if
gotm's avatar
gotm committed
276 277 278 279 280 281
      end do
   end do
   do j=jmin,jmax             ! VEx defined on V-points
      do i=imin,imax
         if (av(i,j) .ge. 1) then
            VEx(i,j)=VEx(i,j)-(PP(i,j+1)-PP(i,j  ))*ARVD1
kbk's avatar
kbk committed
282
         end if
gotm's avatar
gotm committed
283 284
      end do
   end do
285
#endif
gotm's avatar
gotm committed
286 287 288 289 290 291 292 293 294 295 296 297

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

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