uv_diffusion_3d.F90 8.45 KB
Newer Older
gotm's avatar
gotm committed
1 2 3 4
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
hb's avatar
hb committed
5 6
! !ROUTINE: uv_diffusion_3d - hor.\ momentum diffusion
! \label{sec-uv-diffusion-3d}
gotm's avatar
gotm committed
7 8 9 10 11 12
!
! !INTERFACE:
   subroutine uv_diffusion_3d(Am)
!
! !DESCRIPTION:
!
13 14
! Here, the horizontal diffusion terms are discretised in
! momentum-conserving form
hb's avatar
hb committed
15 16
! as well. For simplicity,
! this shown here for these terms as they appear in equations
17 18 19
! (\ref{uEqviCurvi}) and (\ref{vEqviCurvi}),
! i.e.\ without multiplying them by $mn$.
!
hb's avatar
hb committed
20 21 22 23 24 25 26 27 28 29 30 31 32 33
! First horizontal diffusion term in (\ref{uEqviCurvi}):
! \begin{equation}
! \begin{array}{l}
! \displaystyle
! \partial_{\cal X}\left(\frac{2A_Mh_k}{n} m\partial_{\cal X}u_k\right)_{i,j,k}
! \approx
! \\ \\ \displaystyle \quad
! 2A^M_{i+1,j,k} \Delta y_{i+1,j}^c h_{i+1,j,k}^c
! \frac{u_{i+1,j,k}-u_{i,j,k}}{\Delta x_{i+1,j}^c}
! -
! 2A^M_{i,j,k} \Delta y_{i,j}^c h_{i,j,k}^c
! \frac{u_{i,j,k}-u_{i-1,j,k}}{\Delta x_{i,j}^c}
! \end{array}
! \end {equation}
34
!
hb's avatar
hb committed
35 36 37 38 39 40
! Second horizontal diffusion term in (\ref{uEqviCurvi}):
! \begin{equation}
! \begin{array}{l}
! \displaystyle
! \partial_{\cal Y}\left(\frac{A_Mh_k}{m}
! (n\partial_{\cal Y}u_k+m\partial_{\cal X}v_k)\right)_{i,j,k}
41 42
! \approx
! \\ \\ \displaystyle \quad
hb's avatar
hb committed
43 44 45
! \,\,\,\,\,\frac14\left(
! A^M_{i+1,j+1,k}+A^M_{i,j+1,k}+A^M_{i,j,k}+A^M_{i+1,j,k}\right)
! h^+_{i,j,k}\Delta x^+_{i,j}
46
! \\ \\ \displaystyle \quad
hb's avatar
hb committed
47 48 49 50
! \,\,\,\,\,\times\bigg(
! \frac{u_{i,j+1,k}-u_{i,j,k}}{\Delta y^+_{i,j}}
! -\frac{v_{i+1,j,k}-v_{i,j,k}}{\Delta x^+_{i,j}}
! \bigg)
51
! \\ \\ \displaystyle \quad
hb's avatar
hb committed
52 53 54
! -\frac14\left(
! A^M_{i+1,j,k}+A^M_{i,j,k}+A^M_{i,j-1,k}+A^M_{i+1,j-1,k}\right)
! h^+_{i,j-1,k}\Delta x^+_{i,j-1}
55
! \\ \\ \displaystyle \quad
hb's avatar
hb committed
56 57 58 59 60 61
! \,\,\,\,\,\times\bigg(
! \frac{u_{i,j,k}-u_{i,j-1,k}}{\Delta y^+_{i,j-1}}
! -\frac{v_{i+1,j-1,k}-v_{i,j-1,k}}{\Delta x^+_{i,j-1}}
! \bigg)
! \end{array}
! \end {equation}
62
!
hb's avatar
hb committed
63 64 65 66 67 68 69 70 71 72 73 74 75 76
! First horizontal diffusion term in (\ref{vEqviCurvi}):
! \begin{equation}
! \begin{array}{l}
! \displaystyle
! \partial_{\cal Y}\left(\frac{2A_Mh_k}{m} n\partial_{\cal Y}v_k\right)_{i,j,k}
! \approx
! \\ \\ \displaystyle \quad
! 2A^M_{i,j+1,k} \Delta x_{i,j+1}^c h_{i,j+1,k}^c
! \frac{v_{i,j+1,k}-v_{i,j,k}}{\Delta y_{i,j+1}^c}
! -
! 2A^M_{i,j,k} \Delta x_{i,j}^c h_{i,j,k}^c
! \frac{v_{i,j,k}-v_{i,j-1,k}}{\Delta y_{i,j}^c}
! \end{array}
! \end {equation}
77
!
hb's avatar
hb committed
78 79 80 81 82 83
! Second horizontal diffusion term in (\ref{vEqviCurvi}):
! \begin{equation}
! \begin{array}{l}
! \displaystyle
! \partial_{\cal X}\left(\frac{A_Mh_k}{n}
! (n\partial_{\cal Y}u_k+m\partial_{\cal X}v_k)\right)_{i,j,k}
84 85
! \approx
! \\ \\ \displaystyle \quad
hb's avatar
hb committed
86 87 88
! \,\,\,\,\,\frac14\left(
! A^M_{i+1,j+1,k}+A^M_{i,j+1,k}+A^M_{i,j,k}+A^M_{i+1,j,k}\right)
! h^+_{i,j,k}\Delta x^+_{i,j}
89
! \\ \\ \displaystyle \quad
hb's avatar
hb committed
90 91 92 93
! \,\,\,\,\,\times\bigg(
! \frac{u_{i,j+1,k}-u_{i,j,k}}{\Delta y^+_{i,j}}
! -\frac{v_{i+1,j,k}-v_{i,j,k}}{\Delta x^+_{i,j}}
! \bigg)
94
! \\ \\ \displaystyle \quad
hb's avatar
hb committed
95 96 97
! -\frac14\left(
! A^M_{i,j+1,k}+A^M_{i-1,j+1,k}+A^M_{i-1,j,k}+A^M_{i,j,k}\right)
! h^+_{i-1,j,k}\Delta x^+_{i-1,j}
98
! \\ \\ \displaystyle \quad
hb's avatar
hb committed
99 100 101 102 103 104
! \,\,\,\,\,\times\bigg(
! \frac{u_{i-1,j+1,k}-u_{i-1,j,k}}{\Delta y^+_{i-1,j}}
! -\frac{v_{i,j,k}-v_{i-1,j,k}}{\Delta x^+_{i-1,j}}
! \bigg)
! \end{array}
! \end {equation}
105 106
!
! It is assumed here that the horizontal momentum diffusivities
hb's avatar
hb committed
107
! $A^M_{i,j,k}$ are located
108
! on the T-points.
hb's avatar
hb committed
109 110 111 112
!
! For the case of a slice model simulation (compiler option {\tt SLICE\_MODEL}
! activated) the diffusive fluxes in $y$-direction are set to zero.
!
gotm's avatar
gotm committed
113
! !USES:
114
   use domain, only: imin,imax,jmin,jmax,kmax,az,au,av,ax
gotm's avatar
gotm committed
115 116 117 118 119
#if defined(SPHERICAL) || defined(CURVILINEAR)
   use domain, only: dyc,arud1,dxx,dyx,arvd1,dxc
#else
   use domain, only: dx,dy,ard1
#endif
120
   use variables_3d, only: kumin,kvmin,uu,vv,ww,hn,hun,hvn,uuEx,vvEx
bjb's avatar
bjb committed
121
   use getm_timers, only: tic, toc, TIM_UVDIFF3D
bjb's avatar
bjb committed
122
!$ use omp_lib
gotm's avatar
gotm committed
123 124 125 126 127
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
  REALTYPE, intent(in) :: Am

kbk's avatar
kbk committed
128 129 130
! !REVISION HISTORY:
!  Original author(s): Hans Burchard & Karsten Bolding
!
gotm's avatar
gotm committed
131
! !LOCAL VARIABLES:
kbk's avatar
kbk committed
132
   integer                   :: i,j,k,ii,jj,kk
133
   REALTYPE                  :: PP(imin-1:imax+1,jmin-1:jmax+1,1:kmax)
gotm's avatar
gotm committed
134 135 136 137 138 139
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
140
   write(debug,*) 'uv_diffusion_3d() # ',Ncall
gotm's avatar
gotm committed
141
#endif
bjb's avatar
bjb committed
142
   call tic(TIM_UVDIFF3D)
gotm's avatar
gotm committed
143

bjb's avatar
bjb committed
144 145 146
!$OMP PARALLEL DEFAULT(SHARED)                                         &
!$OMP    PRIVATE(i,j,k,ii,jj,kk)

gotm's avatar
gotm committed
147 148
! Central for dx(2*Am*dx(uu^2/hun))
   do k=1,kmax
bjb's avatar
bjb committed
149
!$OMP DO SCHEDULE(RUNTIME)
150 151
      do j=jmin,jmax
         do i=imin,imax+1          ! PP defined on T-points
kbk's avatar
kbk committed
152
            PP(i,j,k)=_ZERO_
gotm's avatar
gotm committed
153 154
            if (az(i,j) .ge. 1) then
               if (k .ge. kumin(i,j)) then
155
                  PP(i,j,k)=_TWO_*Am*DYC*hn(i,j,k)               &
kbk's avatar
kbk committed
156
                      *(uu(i,j,k)/hun(i,j,k)-uu(i-1,j,k)/hun(i-1,j,k))/DXC
gotm's avatar
gotm committed
157 158 159 160
               end if
            end if
         end do
      end do
bjb's avatar
bjb committed
161
!$OMP END DO NOWAIT
gotm's avatar
gotm committed
162
   end do
bjb's avatar
bjb committed
163
!$OMP BARRIER
gotm's avatar
gotm committed
164
   do k=1,kmax
bjb's avatar
bjb committed
165
!$OMP DO SCHEDULE(RUNTIME)
166 167
      do j=jmin,jmax         ! uuEx defined on U-points
         do i=imin,imax
gotm's avatar
gotm committed
168 169 170 171 172 173 174
            if (au(i,j) .ge. 1) then
               if (k .ge. kumin(i,j)) then
                  uuEx(i,j,k)=uuEx(i,j,k)-(PP(i+1,j,k)-PP(i,j,k))*ARUD1
               end if
            end if
         end do
      end do
bjb's avatar
bjb committed
175
!$OMP END DO NOWAIT
gotm's avatar
gotm committed
176
   end do
bjb's avatar
bjb committed
177
!$OMP BARRIER
gotm's avatar
gotm committed
178

179
#ifndef SLICE_MODEL
gotm's avatar
gotm committed
180 181
! Central for dy(Am*(dy(uu^2/hun)+dx(vv^2/hvn)))
   do k=1,kmax
bjb's avatar
bjb committed
182
!$OMP DO SCHEDULE(RUNTIME)
183 184
      do j=jmin-1,jmax          ! PP defined on X-points
         do i=imin,imax
kbk's avatar
kbk committed
185 186
            PP(i,j,k)=_ZERO_
            if (ax(i,j) .ge. 1) then
gotm's avatar
gotm committed
187
               if (k .ge. kumin(i,j)) then
kbk's avatar
kbk committed
188
                  PP(i,j,k)=Am*0.5*(hun(i,j,k)+hun(i,j+1,k))*DXX  &
kbk's avatar
kbk committed
189 190
                      *((uu(i,j+1,k)/hun(i,j+1,k)-uu(i,j,k)/hun(i,j,k))/DYX &
                       +(vv(i+1,j,k)/hvn(i+1,j,k)-vv(i,j,k)/hvn(i,j,k))/DXX )
gotm's avatar
gotm committed
191 192 193 194
               end if
            end if
         end do
      end do
bjb's avatar
bjb committed
195
!$OMP END DO NOWAIT
gotm's avatar
gotm committed
196
   end do
bjb's avatar
bjb committed
197
!$OMP BARRIER
gotm's avatar
gotm committed
198
   do k=1,kmax
bjb's avatar
bjb committed
199
!$OMP DO SCHEDULE(RUNTIME)
200 201
      do j=jmin,jmax
         do i=imin,imax
gotm's avatar
gotm committed
202 203 204 205 206 207 208
            if (au(i,j) .ge. 1) then
               if (k .ge. kumin(i,j)) then
                  uuEx(i,j,k)=uuEx(i,j,k)-(PP(i,j,k)-PP(i,j-1,k))*ARUD1
               end if
            end if
         end do
      end do
bjb's avatar
bjb committed
209
!$OMP END DO NOWAIT
gotm's avatar
gotm committed
210
   end do
bjb's avatar
bjb committed
211
!$OMP BARRIER
212
#endif
gotm's avatar
gotm committed
213 214 215

! Central for dx(Am*(dy(uu^2/hun)+dx(vv^2/hvn)))
   do k=1,kmax
bjb's avatar
bjb committed
216
!$OMP DO SCHEDULE(RUNTIME)
217 218
      do j=jmin,jmax          ! PP defined on X-points
         do i=imin-1,imax
kbk's avatar
kbk committed
219 220
            PP(i,j,k)=_ZERO_
            if (ax(i,j) .ge. 1) then
gotm's avatar
gotm committed
221
               if (k .ge. kumin(i,j)) then
222
                  PP(i,j,k)=Am*_HALF_*(hvn(i+1,j,k)+hvn(i,j,k))*DXX  &
kbk's avatar
kbk committed
223 224
                      *((uu(i,j+1,k)/hun(i,j+1,k)-uu(i,j,k)/hun(i,j,k))/DYX &
                       +(vv(i+1,j,k)/hvn(i+1,j,k)-vv(i,j,k)/hvn(i,j,k))/DXX )
gotm's avatar
gotm committed
225 226 227 228
               end if
            end if
         end do
      end do
bjb's avatar
bjb committed
229
!$OMP END DO NOWAIT
gotm's avatar
gotm committed
230
   end do
bjb's avatar
bjb committed
231
!$OMP BARRIER
gotm's avatar
gotm committed
232
   do k=1,kmax
bjb's avatar
bjb committed
233
!$OMP DO SCHEDULE(RUNTIME)
234 235
      do j=jmin,jmax          ! vvEx defined on V-points
         do i=imin,imax
gotm's avatar
gotm committed
236 237 238 239 240 241 242
            if (av(i,j) .ge. 1) then
               if (k .ge. kvmin(i,j)) then
                  vvEx(i,j,k)=vvEx(i,j,k)-(PP(i,j,k)-PP(i-1,j,k))*ARVD1
               end if
            end if
         end do
      end do
bjb's avatar
bjb committed
243
!$OMP END DO NOWAIT
gotm's avatar
gotm committed
244
   end do
bjb's avatar
bjb committed
245
!$OMP BARRIER
gotm's avatar
gotm committed
246

247
#ifndef SLICE_MODEL
gotm's avatar
gotm committed
248 249
! Central for dy(2*Am*dy(vv^2/hvn))
   do k=1,kmax
bjb's avatar
bjb committed
250
!$OMP DO SCHEDULE(RUNTIME)
251 252
      do j=jmin,jmax+1
         do i=imin,imax          ! PP defined on T-points
gotm's avatar
gotm committed
253 254
            if (az(i,j) .ge. 1) then
               if (k .ge. kvmin(i,j)) then
255
                  PP(i,j,k)=_TWO_*Am*DXC*hn(i,j,k)               &
kbk's avatar
kbk committed
256
                      *(vv(i,j,k)/hvn(i,j,k)-vv(i,j-1,k)/hvn(i,j-1,k))/DYC
gotm's avatar
gotm committed
257 258 259 260
               end if
            end if
         end do
      end do
bjb's avatar
bjb committed
261
!$OMP END DO NOWAIT
gotm's avatar
gotm committed
262
   end do
bjb's avatar
bjb committed
263
!$OMP BARRIER
gotm's avatar
gotm committed
264 265

   do k=1,kmax
bjb's avatar
bjb committed
266
!$OMP DO SCHEDULE(RUNTIME)
267 268
      do j=jmin,jmax          ! vvEx defined on V-points
         do i=imin,imax
gotm's avatar
gotm committed
269 270 271 272 273 274 275
            if (av(i,j) .ge. 1) then
               if (k .ge. kvmin(i,j)) then
                  vvEx(i,j,k)=(vvEx(i,j,k)-(PP(i,j+1,k)-PP(i,j,k))*ARVD1)
               end if
            end if
         end do
      end do
bjb's avatar
bjb committed
276
!$OMP END DO NOWAIT
gotm's avatar
gotm committed
277
   end do
278
#endif
bjb's avatar
bjb committed
279
!$OMP END PARALLEL
gotm's avatar
gotm committed
280

bjb's avatar
bjb committed
281
   call toc(TIM_UVDIFF3D)
gotm's avatar
gotm committed
282 283 284 285 286 287 288 289 290 291 292
#ifdef DEBUG
   write(debug,*) 'Leaving uv_diffusion_3d()'
   write(debug,*)
#endif
   return
   end subroutine uv_diffusion_3d
!EOC

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