uv_advect.F90 7.42 KB
Newer Older
kbk's avatar
kbk committed
1
!$Id: uv_advect.F90,v 1.10 2006-03-01 15:54:07 kbk Exp $
gotm's avatar
gotm committed
2 3 4 5
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
6
! !IROUTINE: uv_advect - 2D advection of momentum \label{sec-uv-advect}
gotm's avatar
gotm committed
7 8 9 10 11 12
!
! !INTERFACE:
   subroutine uv_advect
!
! !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 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
! The advective terms in the vertically integrated 
! momentum equation are discretised in
! a momentum-conservative form. This is carried out here for the 
! advective terms in the $U$-equation (\ref{UMOM}) and the 
! $V$-equation (\ref{VMOM}) (after applying the curvilinear
! coordinate transformationand multiplying these
! equations with $mn$). 
! 
! First advection term in (\ref{UMOM}):
! \begin{equation}
! \begin{array}{l}
! \displaystyle
! \left(mn\,\partial_{\cal X}\left(\frac{U^2}{Dn}\right)\right)_{i,j}\approx \\ \\
! \quad
! \displaystyle
! \frac{
! \frac12(U_{i+1,j}+U_{i,j})\tilde u_{i+1,j}\Delta y^c_{i+1,j}-
! \frac12(U_{i,j}+U_{i-1,j})\tilde u_{i,j}\Delta y^c_{i,j}
! }{\Delta x^u_{i,j}\Delta y^u_{i,j}}
! \end{array}
! \end{equation}
! 
! For the upwind scheme used here, the inter-facial velocities which are defined
! on T-points are here
! calculated as:
! 
! \begin{equation}
! \tilde u_{i,j}=
! \left\{
! \begin{array}{ll}
! \displaystyle
! \frac{U_{i-1,j}}{D^u_{i-1,j}} & \mbox{ for } \frac12(U_{i,j}+U_{i-1,j})>0\\ \\
! \displaystyle
! \frac{U_{i,j}}{D^u_{i,j}} & \mbox{ else. } 
! \end{array}
! \right.
! \end{equation}
! 
! Second advection term in (\ref{UMOM}):
! \begin{equation}
! \begin{array}{l}
! \displaystyle 
! \left(mn\,\partial_{\cal Y}y\left(\frac{UV}{Dm}\right)\right)_{i,j,k}\approx \\ \\ 
! \displaystyle 
! \quad
! \frac{
! \frac12(V_{i+1,j}+V_{i,j})\tilde u_{i,j}\Delta x^+_{i,j}-
! \frac12(V_{i+1,j-1}+V_{i,j-1})\tilde u_{i,j-1}\Delta x^+_{i,j-1}
! }{\Delta x^u_{i,j}\Delta y^u_{i,j}}
! \end{array}
! \end{equation}
! 
! For the upwind scheme used here, the inter-facial 
! velocities which are defined on
! X-points are here
! calculated as:
! 
! \begin{equation}
! \tilde u_{i,j}=
! \left\{
! \begin{array}{ll}
! \displaystyle
! \frac{U_{i,j}}{D^u_{i,j}} & \mbox{ for } \frac12(V_{i+1,j,k}+V_{i,j,k})>0\\ \\
! \displaystyle
! \frac{U_{i,j+1}}{D^u_{i,j+1}} & \mbox{ else. } 
! \end{array}
! \right.
! \end{equation}
! 
! First advection term in (\ref{VMOM}):
! \begin{equation}
! \begin{array}{l}
! \displaystyle 
! \left(mn\,\partial_{\cal X}\left(\frac{UV}{Dn}\right)\right)_{i,j,k}\approx \\ \\ 
! \displaystyle 
! \quad
! \frac{
! \frac12(U_{i,j+1}+U_{i,j})\tilde v_{i,j}\Delta y^+_{i,j}-
! \frac12(U_{i-1,j+1}+U_{i-1,j})\tilde v_{i-1,j}\Delta y^+_{i-1,j}
! }{\Delta x^v_{i,j}\Delta y^v_{i,j}}
! \end{array}
! \end{equation}
! 
! For the upwind scheme used here, the interfacial 
! velocities which are defined on
! X-points are here
! calculated as:
! 
! \begin{equation}
! \tilde v_{i,j}=
! \left\{
! \begin{array}{ll}
! \displaystyle
! \frac{V_{i,j}}{D^v_{i,j}} & \mbox{ for } \frac12(U_{i+1,j}+U_{i,j})>0\\ \\
! \displaystyle
! \frac{V_{i+1,j}}{D^v_{i+1,j}} & \mbox{ else. } 
! \end{array}
! \right.
! \end{equation}
! 
! Second advection term in (\ref{VMOM}):
! \begin{equation}
! \begin{array}{l}
! \displaystyle
! \left(mn\,\partial_{\cal Y}\left(\frac{V^2}{Dm}\right)\right)_{i,j,k}\approx \\ \\
! \quad
! \displaystyle
! \frac{
! \frac12(V_{i,j+1}+V_{i,j})\tilde v_{i,j+1}\Delta x^c_{i,j+1}-
! \frac12(V_{i,j}+V_{i,j-1})\tilde v_{i,j}\Delta x^c_{i,j}
! }{\Delta x^v_{i,j}\Delta y^v_{i,j}}
! \end{array}
! \end{equation}
! 
! For the upwind scheme used here, the interfacial velocities which are defined
! on T-points are here
! calculated as:
! 
! \begin{equation}
! \tilde v_{i,j}=
! \left\{
! \begin{array}{ll}
! \displaystyle
! \frac{V_{i,j-1}}{D^v_{i,j-1}} & \mbox{ for } \frac12(V_{i,j}+V_{i,j-1})>0\\ \\
! \displaystyle
! \frac{V_{i,j}}{D^v_{i,j}} & \mbox{ else. } 
! \end{array}
! \right.
! \end{equation}
!
! When working with the option {\tt SLICE\_MODEL}, the calculation of
! all gradients in $y$-direction is suppressed.
!
gotm's avatar
gotm committed
146
! !USES:
kbk's avatar
kbk committed
147
   use domain, only: imin,imax,jmin,jmax,az,au,av,ax
kbk's avatar
kbk committed
148
   use domain, only: ioff,joff
gotm's avatar
gotm committed
149 150 151 152 153
#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
154
   use variables_2d, only: U,DU,UEx,V,DV,VEx,PP
gotm's avatar
gotm committed
155 156 157 158 159 160 161 162
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
kbk's avatar
kbk committed
163 164 165
! !REVISION HISTORY:
!  Original author(s): Hans Burchard & Karsten Bolding
!
gotm's avatar
gotm committed
166
! !LOCAL VARIABLES:
kbk's avatar
kbk committed
167
   integer                   :: i,j,ii,jj
gotm's avatar
gotm committed
168 169 170 171 172 173 174 175 176
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'uv_advect() # ',Ncall
#endif

kbk's avatar
kbk committed
177 178
!  Upstream for dx(U^2/D)
   do j=jmin,jmax         ! PP defined on T-points
gotm's avatar
gotm committed
179
      do i=imin,imax+1
kbk's avatar
kbk committed
180
         if (az(i,j) .ge. 1) then
gotm's avatar
gotm committed
181
            PP(i,j)=0.5*(U(i-1,j)+U(i,j))
kbk's avatar
kbk committed
182
            if (PP(i,j) .gt. _ZERO_) then
kbk's avatar
kbk committed
183
               ii=i-1
kbk's avatar
kbk committed
184
            else
kbk's avatar
kbk committed
185 186 187
               ii=i
            end if
            PP(i,j)=PP(i,j)*U(ii,j)/DU(ii,j)*DYC
kbk's avatar
kbk committed
188 189
         else
            PP(i,j) = _ZERO_
kbk's avatar
kbk committed
190
         end if
gotm's avatar
gotm committed
191 192 193 194
      end do
   end do
   do j=jmin,jmax      ! UEx defined on U-points
      do i=imin,imax
hb's avatar
hb committed
195
         if (au(i,j) .eq. 1) then
gotm's avatar
gotm committed
196
            UEx(i,j)=(PP(i+1,j)-PP(i  ,j))*ARUD1
kbk's avatar
kbk committed
197
         end if
gotm's avatar
gotm committed
198 199 200
      end do
   end do

201
#ifndef SLICE_MODEL
kbk's avatar
kbk committed
202
!  Upstream for dy(UV/D)
gotm's avatar
gotm committed
203
   do j=jmin-1,jmax        ! PP defined on X-points
kbk's avatar
kbk committed
204 205
      do i=imin,imax
         if (ax(i,j) .ge. 1) then
gotm's avatar
gotm committed
206
            PP(i,j)=0.5*(V(i+1,j)+V(i,j))
kbk's avatar
kbk committed
207 208 209 210 211 212
            if (PP(i,j) .gt. _ZERO_) then
               jj=j
            else
               jj=j+1
            end if
            PP(i,j)=PP(i,j)*U(i,jj)/DU(i,jj)*DXX
kbk's avatar
kbk committed
213 214
         else
            PP(i,j) = _ZERO_
kbk's avatar
kbk committed
215
         end if
gotm's avatar
gotm committed
216 217 218 219
      end do
   end do
   do j=jmin,jmax        !UEx defined on U-points
      do i=imin,imax
hb's avatar
hb committed
220
         if (au(i,j) .eq. 1) then
gotm's avatar
gotm committed
221
            UEx(i,j)=UEx(i,j)+(PP(i,j  )-PP(i,j-1))*ARUD1
kbk's avatar
kbk committed
222
         end if
gotm's avatar
gotm committed
223 224
      end do
   end do
225
#endif
gotm's avatar
gotm committed
226 227

! Upstream for dx(UV/D)
kbk's avatar
kbk committed
228
   do j=jmin,jmax      ! PP defined on X-points
gotm's avatar
gotm committed
229
      do i=imin-1,imax
kbk's avatar
kbk committed
230
         if (ax(i,j) .ge. 1) then
gotm's avatar
gotm committed
231
            PP(i,j)=0.5*(U(i,j)+U(i,j+1))
kbk's avatar
kbk committed
232 233 234 235 236 237
            if (PP(i,j) .gt. _ZERO_) then
               ii=i
            else
               ii=i+1
            end if
            PP(i,j)=PP(i,j)*V(ii,j)/DV(ii,j)*DYX
kbk's avatar
kbk committed
238 239
         else
            PP(i,j) = _ZERO_
kbk's avatar
kbk committed
240
         end if
gotm's avatar
gotm committed
241 242 243 244
      end do
   end do
   do j=jmin,jmax          ! VEx defined on V-points
      do i=imin,imax
hb's avatar
hb committed
245
         if (av(i,j) .eq. 1) then
gotm's avatar
gotm committed
246
            VEx(i,j)=(PP(i  ,j)-PP(i-1,j))*ARVD1
kbk's avatar
kbk committed
247
         end if
gotm's avatar
gotm committed
248 249 250
      end do
   end do

251
#ifndef SLICE_MODEL
kbk's avatar
kbk committed
252
!  Upstream for dy(V^2/D)
gotm's avatar
gotm committed
253
   do j=jmin,jmax+1     ! PP defined on T-points
kbk's avatar
kbk committed
254 255
      do i=imin,imax
         if (az(i,j) .ge. 1) then
gotm's avatar
gotm committed
256
            PP(i,j)=0.5*(V(i,j-1)+V(i,j))
kbk's avatar
kbk committed
257 258 259 260 261 262
            if (PP(i,j) .gt. _ZERO_) then
               jj=j-1
            else
               jj=j
            end if
            PP(i,j)=PP(i,j)*V(i,jj)/DV(i,jj)*DXC
kbk's avatar
kbk committed
263 264
         else
            PP(i,j) = _ZERO_
kbk's avatar
kbk committed
265
         end if
gotm's avatar
gotm committed
266 267 268 269
      end do
   end do
   do j=jmin,jmax             ! VEx defined on V-points
      do i=imin,imax
hb's avatar
hb committed
270
         if (av(i,j) .eq. 1) then
gotm's avatar
gotm committed
271
            VEx(i,j)=VEx(i,j)+(PP(i,j+1)-PP(i,j  ))*ARVD1
kbk's avatar
kbk committed
272
         end if
gotm's avatar
gotm committed
273 274
      end do
   end do
275
#endif
gotm's avatar
gotm committed
276 277 278 279 280 281 282 283 284 285 286 287

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

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