uv_advect.F90 7.78 KB
Newer Older
gotm's avatar
gotm committed
1 2 3 4
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
5
! !IROUTINE: uv_advect - 2D advection of momentum \label{sec-uv-advect}
gotm's avatar
gotm committed
6 7 8 9 10 11
!
! !INTERFACE:
   subroutine uv_advect
!
! !DESCRIPTION:
!
hb's avatar
hb committed
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
! 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
145
! !USES:
kbk's avatar
kbk committed
146
   use domain, only: imin,imax,jmin,jmax,az,au,av,ax
kbk's avatar
kbk committed
147
   use domain, only: ioff,joff
gotm's avatar
gotm committed
148 149 150 151 152
#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
153
   use variables_2d, only: U,DU,UEx,V,DV,VEx,PP
bjb's avatar
bjb committed
154
   use getm_timers, only: tic, toc, TIM_UVADVECT
bjb's avatar
bjb committed
155
!$ use omp_lib
gotm's avatar
gotm committed
156 157
   IMPLICIT NONE
!
kbk's avatar
kbk committed
158 159 160
! !REVISION HISTORY:
!  Original author(s): Hans Burchard & Karsten Bolding
!
gotm's avatar
gotm committed
161
! !LOCAL VARIABLES:
kbk's avatar
kbk committed
162
   integer                   :: i,j,ii,jj
gotm's avatar
gotm committed
163 164 165 166 167 168 169 170
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'uv_advect() # ',Ncall
#endif
bjb's avatar
bjb committed
171
   CALL tic(TIM_UVADVECT)
gotm's avatar
gotm committed
172

bjb's avatar
bjb committed
173 174
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,ii,jj)

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

203
#ifndef SLICE_MODEL
kbk's avatar
kbk committed
204
!  Upstream for dy(UV/D)
bjb's avatar
bjb committed
205
!$OMP DO SCHEDULE(RUNTIME)
gotm's avatar
gotm committed
206
   do j=jmin-1,jmax        ! PP defined on X-points
kbk's avatar
kbk committed
207 208
      do i=imin,imax
         if (ax(i,j) .ge. 1) then
gotm's avatar
gotm committed
209
            PP(i,j)=0.5*(V(i+1,j)+V(i,j))
kbk's avatar
kbk committed
210 211 212 213 214 215
            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
216 217
         else
            PP(i,j) = _ZERO_
kbk's avatar
kbk committed
218
         end if
gotm's avatar
gotm committed
219 220
      end do
   end do
bjb's avatar
bjb committed
221 222
!$OMP END DO
!$OMP DO SCHEDULE(RUNTIME)
gotm's avatar
gotm committed
223 224
   do j=jmin,jmax        !UEx defined on U-points
      do i=imin,imax
hb's avatar
hb committed
225
         if (au(i,j) .eq. 1) then
gotm's avatar
gotm committed
226
            UEx(i,j)=UEx(i,j)+(PP(i,j  )-PP(i,j-1))*ARUD1
kbk's avatar
kbk committed
227
         end if
gotm's avatar
gotm committed
228 229
      end do
   end do
bjb's avatar
bjb committed
230
!$OMP END DO
231
#endif
gotm's avatar
gotm committed
232 233

! Upstream for dx(UV/D)
bjb's avatar
bjb committed
234
!$OMP DO SCHEDULE(RUNTIME)
kbk's avatar
kbk committed
235
   do j=jmin,jmax      ! PP defined on X-points
gotm's avatar
gotm committed
236
      do i=imin-1,imax
kbk's avatar
kbk committed
237
         if (ax(i,j) .ge. 1) then
gotm's avatar
gotm committed
238
            PP(i,j)=0.5*(U(i,j)+U(i,j+1))
kbk's avatar
kbk committed
239 240 241 242 243 244
            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
245 246
         else
            PP(i,j) = _ZERO_
kbk's avatar
kbk committed
247
         end if
gotm's avatar
gotm committed
248 249
      end do
   end do
bjb's avatar
bjb committed
250 251
!$OMP END DO
!$OMP DO SCHEDULE(RUNTIME)
gotm's avatar
gotm committed
252 253
   do j=jmin,jmax          ! VEx defined on V-points
      do i=imin,imax
hb's avatar
hb committed
254
         if (av(i,j) .eq. 1) then
gotm's avatar
gotm committed
255
            VEx(i,j)=(PP(i  ,j)-PP(i-1,j))*ARVD1
kbk's avatar
kbk committed
256
         end if
gotm's avatar
gotm committed
257 258
      end do
   end do
bjb's avatar
bjb committed
259
!$OMP END DO
gotm's avatar
gotm committed
260

261
#ifndef SLICE_MODEL
kbk's avatar
kbk committed
262
!  Upstream for dy(V^2/D)
bjb's avatar
bjb committed
263
!$OMP DO SCHEDULE(RUNTIME)
gotm's avatar
gotm committed
264
   do j=jmin,jmax+1     ! PP defined on T-points
kbk's avatar
kbk committed
265 266
      do i=imin,imax
         if (az(i,j) .ge. 1) then
gotm's avatar
gotm committed
267
            PP(i,j)=0.5*(V(i,j-1)+V(i,j))
kbk's avatar
kbk committed
268 269 270 271 272 273
            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
274 275
         else
            PP(i,j) = _ZERO_
kbk's avatar
kbk committed
276
         end if
gotm's avatar
gotm committed
277 278
      end do
   end do
bjb's avatar
bjb committed
279 280
!$OMP END DO
!$OMP DO SCHEDULE(RUNTIME)
gotm's avatar
gotm committed
281 282
   do j=jmin,jmax             ! VEx defined on V-points
      do i=imin,imax
hb's avatar
hb committed
283
         if (av(i,j) .eq. 1) then
gotm's avatar
gotm committed
284
            VEx(i,j)=VEx(i,j)+(PP(i,j+1)-PP(i,j  ))*ARVD1
kbk's avatar
kbk committed
285
         end if
gotm's avatar
gotm committed
286 287
      end do
   end do
bjb's avatar
bjb committed
288
!$OMP END DO
289
#endif
gotm's avatar
gotm committed
290

bjb's avatar
bjb committed
291 292
!$OMP END PARALLEL

bjb's avatar
bjb committed
293
   CALL toc(TIM_UVADVECT)
gotm's avatar
gotm committed
294 295 296 297 298 299 300 301 302 303 304
#ifdef DEBUG
     write(debug,*) 'Leaving uv_advect()'
     write(debug,*)
#endif
   return
   end subroutine uv_advect
!EOC

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