vv_momentum_3d.F90 14.9 KB
Newer Older
gotm's avatar
gotm committed
1 2 3 4
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
hb's avatar
hb committed
5
! !ROUTINE: vv_momentum_3d - $y$-momentum eq.\ \label{sec-vv-momentum-3d}
gotm's avatar
gotm committed
6 7
!
! !INTERFACE:
8
   subroutine vv_momentum_3d(n,bdy3d)
gotm's avatar
gotm committed
9 10 11
!
! !DESCRIPTION:
!
hb's avatar
hb committed
12 13 14 15 16 17 18 19 20
! Here, the budget equation for layer-averaged momentum in eastern direction,
! $q_k$,
! is calculated. The physical equation is given as equation (\ref{vEq}),
! the layer-integrated equation as (\ref{vEqvi}), and after curvilinear
! transformation as (\ref{vEqviCurvi}).
! In this routine, first the Coriolis rotation term, $fp_k$ is calculated,
! either as direct transport averaging, or following \cite{ESPELIDea00}
! by using velocity averages (in case the compiler option {\tt NEW\_CORI}
! is set).
21
!
hb's avatar
hb committed
22 23 24 25 26 27 28 29
! As a next step, explicit forcing terms (advection, diffusion,
! internal pressure gradient, surface stresses) are added up (into the variable
! {\tt ex(k)}), the eddy viscosity is horizontally interpolated to the V-point,
! and the barotropic pressure gradient is calculated (the latter
! includes the pressure gradient correction for drying points, see
! section \ref{Section_dry}).
! Afterwards, the matrix is set up for each water column, and it is solved
! by means of a tri-diagonal matrix solver.
30 31 32 33 34 35 36 37 38
!
! In case that the compiler option {\tt STRUCTURE\_FRICTION} is switched on,
! the frictional effect of structures in the water column is calculated
! by adding the quadratic frictional term $C v \sqrt{u^2+v^2}$ (with a minus
! sign on
! the right hand side) numerically implicitly to the $v$-equation,
! with the friction coefficient $C$. The explicit part of this term,
! $C\sqrt{u^2+v^2}$,
! is calculated in the routine {\tt structure\_friction\_3d.F90}.
39
!
hb's avatar
hb committed
40 41 42
! Finally, the new velocity profile is shifted such that its vertical
! integral is identical to the time integral of the vertically integrated
! transport.
43
!
hb's avatar
hb committed
44 45 46 47 48
! When GETM is run as a slice model (compiler option {\tt SLICE\_MODEL}
! is activated), the result for $j=2$ is copied to $j=1$ and $j=3$.
! If the compiler option {\tt XZ\_PLUME\_TEST} is set, a slope
! of {\tt yslope} for bottom and isopycnals into the $y$-direction is
! prescribed, which has to be hard-coded as local variable.
49
!
gotm's avatar
gotm committed
50
! !USES:
51
   use exceptions
gotm's avatar
gotm committed
52
   use parameters, only: g,avmmol,rho_0
Knut's avatar
Knut committed
53
   use domain, only: imin,imax,jmin,jmax,kmax,H,min_depth
kbk's avatar
kbk committed
54
   use domain, only: dry_v,corv,au,av,az,ax
gotm's avatar
gotm committed
55
#if defined CURVILINEAR || defined SPHERICAL
56
   use domain, only: dyv,arvd1,dxc,dyx,dyc,dxx
gotm's avatar
gotm committed
57
#else
58
   use domain, only: dx,dy
gotm's avatar
gotm committed
59
#endif
60 61
   use variables_2d, only: Vint,D
   use bdy_3d, only: do_bdy_3d
gotm's avatar
gotm committed
62
   use variables_3d, only: dt,cnpar,kvmin,uu,vv,huo,hvo,hvn,vvEx,ww,hun
Knut's avatar
Knut committed
63
   use variables_3d, only: num,nuh,sseo,Dvn,rrv
Hans Burchard's avatar
Hans Burchard committed
64
#ifdef _MOMENTUM_TERMS_
65
   use variables_3d, only: tdv_v,cor_v,ipg_v,epg_v,vsd_v,hsd_v
Hans Burchard's avatar
Hans Burchard committed
66
#endif
67
#ifdef XZ_PLUME_TEST
68
   use variables_3d, only: buoy
69
#endif
70 71 72
#ifdef STRUCTURE_FRICTION
   use variables_3d, only: sf
#endif
73 74
#ifndef NO_BAROCLINIC
   use variables_3d, only: idpdy
gotm's avatar
gotm committed
75
#endif
76 77
   use halo_zones, only: update_3d_halo,wait_halo,V_TAG
   use meteo, only: tausy,airp
78
   use m3d, only: ip_fac
79
   use m3d, only: vel_check,min_vel,max_vel
bjb's avatar
bjb committed
80
   use getm_timers, only: tic, toc, TIM_VVMOMENTUM, TIM_VVMOMENTUMH
bjb's avatar
bjb committed
81
!$ use omp_lib
gotm's avatar
gotm committed
82 83 84
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
85
   integer, intent(in)                 :: n
kbk's avatar
kbk committed
86
   logical, intent(in)                 :: bdy3d
gotm's avatar
gotm committed
87
!
kbk's avatar
kbk committed
88 89 90
! !REVISION HISTORY:
!  Original author(s): Hans Burchard & Karsten Bolding
!
gotm's avatar
gotm committed
91
! !LOCAL VARIABLES:
kbk's avatar
kbk committed
92
   integer                   :: i,j,k,rc
bjb's avatar
bjb committed
93 94 95 96 97
   REALTYPE, POINTER         :: dif(:)
   REALTYPE, POINTER         :: auxn(:),auxo(:)
   REALTYPE, POINTER         :: a1(:),a2(:)
   REALTYPE, POINTER         :: a3(:),a4(:)
   REALTYPE, POINTER         :: Res(:),ex(:)
kbk's avatar
kbk committed
98 99 100
   REALTYPE                  :: zp,zm,zy,ResInt,Diff,Uloc
   REALTYPE                  :: gamma=g*rho_0
   REALTYPE                  :: cord_curv=_ZERO_
bjb's avatar
bjb committed
101
   REALTYPE                  :: gammai,rho_0i
hb's avatar
hb committed
102 103 104
#ifdef XZ_PLUME_TEST
   REALTYPE                  :: yslope=0.001
#endif
105
   integer                   :: status
gotm's avatar
gotm committed
106 107 108 109 110 111 112 113
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'vv_momentum_3d() # ',Ncall
#endif
bjb's avatar
bjb committed
114 115
   call tic(TIM_VVMOMENTUM)

bjb's avatar
bjb committed
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
   gammai=_ONE_/gamma
   rho_0i=_ONE_/rho_0


!$OMP PARALLEL DEFAULT(SHARED)                                         &
!$OMP    PRIVATE(i,j,k,rc,zp,zm,zy,ResInt,Diff,Uloc,cord_curv)         &
!$OMP    PRIVATE(dif,auxn,auxo,a1,a2,a3,a4,Res,ex)

 ! Each thread allocates its own HEAP storage:
   allocate(dif(1:kmax-1),stat=rc)    ! work array
   if (rc /= 0) stop 'vv_momentum_3d: Error allocating memory (dif)'

   allocate(auxn(1:kmax-1),stat=rc)    ! work array
   if (rc /= 0) stop 'vv_momentum_3d: Error allocating memory (auxn)'

   allocate(auxo(1:kmax-1),stat=rc)    ! work array
   if (rc /= 0) stop 'vv_momentum_3d: Error allocating memory (auxo)'

   allocate(a1(0:kmax),stat=rc)    ! work array
   if (rc /= 0) stop 'vv_momentum_3d: Error allocating memory (a1)'

   allocate(a2(0:kmax),stat=rc)    ! work array
   if (rc /= 0) stop 'vv_momentum_3d: Error allocating memory (a2)'

   allocate(a3(0:kmax),stat=rc)    ! work array
   if (rc /= 0) stop 'vv_momentum_3d: Error allocating memory (a3)'

   allocate(a4(0:kmax),stat=rc)    ! work array
   if (rc /= 0) stop 'vv_momentum_3d: Error allocating memory (a4)'

   allocate(Res(0:kmax),stat=rc)    ! work array
   if (rc /= 0) stop 'vv_momentum_3d: Error allocating memory (Res)'

   allocate(ex(0:kmax),stat=rc)    ! work array
   if (rc /= 0) stop 'vv_momentum_3d: Error allocating memory (ex)'

! Note: We do not need to initialize these work arrays.
153
!   Tested BJB 2009-09-25.
bjb's avatar
bjb committed
154 155

!$OMP DO SCHEDULE(RUNTIME)
156 157
   do j=jmin,jmax
      do i=imin,imax
gotm's avatar
gotm committed
158 159

         if ((av(i,j) .eq. 1) .or. (av(i,j) .eq. 2)) then
kbk's avatar
kbk committed
160

gotm's avatar
gotm committed
161 162 163 164 165 166 167 168 169
            if (kmax .gt. kvmin(i,j)) then

               do k=kvmin(i,j),kmax      ! explicit terms
! Espelid et al. [2000], IJNME 49, 1521-1545
#ifdef NEW_CORI
                Uloc=(uu(i  ,j  ,k)/sqrt(huo(i  ,j  ,k))  &
                     +uu(i-1,j  ,k)/sqrt(huo(i-1,j  ,k))  &
                     +uu(i  ,j+1,k)/sqrt(huo(i  ,j+1,k))  &
                     +uu(i-1,j+1,k)/sqrt(huo(i-1,j+1,k))) &
bjb's avatar
bjb committed
170
                     *_QUART_*sqrt(hvo(i,j,k))
gotm's avatar
gotm committed
171
#else
bjb's avatar
bjb committed
172
                Uloc=_QUART_*(uu(i,j,k)+uu(i-1,j,k)+uu(i,j+1,k)+uu(i-1,j+1,k))
gotm's avatar
gotm committed
173
#endif
174 175
#if defined(SPHERICAL) || defined(CURVILINEAR)
                  cord_curv=(vv(i,j,k)*(DYX-DYXIM1)-Uloc*(DXCJP1-DXC))     &
gotm's avatar
gotm committed
176
                        /hvo(i,j,k)*ARVD1
kbk's avatar
kbk committed
177
                  ex(k)=(cord_curv-corv(i,j))*Uloc
178
#else
kbk's avatar
kbk committed
179
                  ex(k)=-corv(i,j)*Uloc
180
#endif
Hans Burchard's avatar
Hans Burchard committed
181 182 183
#ifdef _MOMENTUM_TERMS_
                  cor_v(i,j,k)=-dry_v(i,j)*ex(k)
#endif
184 185
#ifdef NO_BAROCLINIC
                  ex(k)=dry_v(i,j)*(ex(k)-vvEx(i,j,k))
186 187
#else
#ifdef XZ_PLUME_TEST
188
                  ex(k)=dry_v(i,j)*(ex(k)-vvEx(i,j,k)+idpdy(i,j,k)+yslope*hvn(i,j,k)*(buoy(i,j,kmax)-buoy(i,j,k)))
gotm's avatar
gotm committed
189
#else
190
                  ex(k)=dry_v(i,j)*(ex(k)-vvEx(i,j,k)+ip_fac*idpdy(i,j,k))
191
#endif
Hans Burchard's avatar
Hans Burchard committed
192 193 194
#ifdef _MOMENTUM_TERMS_
                  ipg_v(i,j,k)=-dry_v(i,j)*ip_fac*idpdy(i,j,k)
#endif
gotm's avatar
gotm committed
195 196 197
#endif
               end do
               ex(kmax)=ex(kmax)                                      &
Hans Burchard's avatar
Hans Burchard committed
198
                       +dry_v(i,j)*_HALF_*(tausy(i,j)+tausy(i,j+1))*rho_0i
gotm's avatar
gotm committed
199 200
!     Eddy viscosity
               do k=kvmin(i,j),kmax-1
bjb's avatar
bjb committed
201
                  dif(k)=_HALF_*(num(i,j,k)+num(i,j+1,k)) + avmmol
gotm's avatar
gotm committed
202 203 204 205 206
               end do

!     Auxiliury terms, old and new time level,
!     cnpar: Crank-Nicholson parameter
               do k=kvmin(i,j),kmax-1
207 208
                  auxo(k)=2*(1-cnpar)*dt*dif(k)/(hvo(i,j,k+1)+hvo(i,j,k))
                  auxn(k)=2*   cnpar *dt*dif(k)/(hvn(i,j,k+1)+hvn(i,j,k))
gotm's avatar
gotm committed
209 210 211
               end do

!     Barotropic pressure gradient
212
#ifndef NO_BAROTROPIC
gotm's avatar
gotm committed
213 214 215
               zp=max(sseo(i,j+1),-H(i,j  )+min(min_depth,D(i,j+1)))
               zm=max(sseo(i,j  ),-H(i,j+1)+min(min_depth,D(i,j  )))
               zy=(zp-zm+(airp(i,j+1)-airp(i,j))/gamma)/DYV
216
#else
bjb's avatar
bjb committed
217
               zy=_ZERO_
218
#endif
gotm's avatar
gotm committed
219 220 221 222 223

!     Matrix elements for surface layer
               k=kmax
               a1(k)=-auxn(k-1)/hvn(i,j,k-1)
               a2(k)=1+auxn(k-1)/hvn(i,j,k)
kbk's avatar
kbk committed
224 225 226
               a4(k)=vv(i,j,k  )*(1-auxo(k-1)/hvo(i,j,k))              &
                    +vv(i,j,k-1)*auxo(k-1)/hvo(i,j,k-1)                &
                    +dt*ex(k)                                          &
bjb's avatar
bjb committed
227
                    -dt*_HALF_*(hvo(i,j,k)+hvn(i,j,k))*g*zy
gotm's avatar
gotm committed
228 229 230 231 232

!     Matrix elements for inner layers
               do k=kvmin(i,j)+1,kmax-1
                  a3(k)=-auxn(k  )/hvn(i,j,k+1)
                  a1(k)=-auxn(k-1)/hvn(i,j,k-1)
bjb's avatar
bjb committed
233
                  a2(k)=_ONE_+(auxn(k)+auxn(k-1))/hvn(i,j,k)
kbk's avatar
kbk committed
234 235 236 237
                  a4(k)=vv(i,j,k+1)*auxo(k)/hvo(i,j,k+1)               &
                       +vv(i,j,k  )*(1-(auxo(k)+auxo(k-1))/hvo(i,j,k)) &
                       +vv(i,j,k-1)*auxo(k-1)/hvo(i,j,k-1)             &
                       +dt*ex(k)                                       &
bjb's avatar
bjb committed
238
                       -dt*_HALF_*(hvo(i,j,k)+hvn(i,j,k))*g*zy
gotm's avatar
gotm committed
239 240 241 242 243
               end do

!     Matrix elements for bottom layer
               k=kvmin(i,j)
               a3(k)=-auxn(k  )/hvn(i,j,k+1)
bjb's avatar
bjb committed
244 245
               a2(k)=_ONE_+auxn(k)/hvn(i,j,k)                          &
                        +dt*rrv(i,j)/(_HALF_*(hvn(i,j,k)+hvo(i,j,k)))
kbk's avatar
kbk committed
246
               a4(k)=vv(i,j,k+1)*auxo(k)/hvo(i,j,k+1)                  &
bjb's avatar
bjb committed
247
                       +vv(i,j,k  )*(_ONE_-auxo(k)/hvo(i,j,k))         &
kbk's avatar
kbk committed
248
                       +dt*ex(k)                                       &
bjb's avatar
bjb committed
249
                       -dt*_HALF_*(hvo(i,j,k)+hvn(i,j,k))*g*zy
gotm's avatar
gotm committed
250

251 252
#ifdef STRUCTURE_FRICTION
               do k=kvmin(i,j),kmax
bjb's avatar
bjb committed
253
                  a2(k)=a2(k)+dt*_HALF_*(sf(i,j,k)+sf(i,j+1,k))
254 255 256
               end do
#endif

gotm's avatar
gotm committed
257 258 259 260 261 262 263 264 265
               call getm_tridiagonal(kmax,kvmin(i,j),kmax,a1,a2,a3,a4,Res)

!     Transport correction: the integral of the new velocities has to
!     be the same than the transport calculated by the external mode, Vint.

               ResInt= _ZERO_
               do k=kvmin(i,j),kmax
                  ResInt=ResInt+Res(k)
               end do
Knut's avatar
Knut committed
266
               Diff=(Vint(i,j)-ResInt)/Dvn(i,j)
gotm's avatar
gotm committed
267 268

               do k=kvmin(i,j),kmax
Hans Burchard's avatar
Hans Burchard committed
269 270 271 272
#ifdef _MOMENTUM_TERMS_
                  tdv_v(i,j,k)=vv(i,j,k)
                  epg_v(i,j,k)=_HALF_*(hvo(i,j,k)+hvn(i,j,k))*g*zy     &
                              -hvn(i,j,k)*Diff/dt
273
                  hsd_v(i,j,k)=hsd_v(i,j,k)*dry_v(i,j)
Hans Burchard's avatar
Hans Burchard committed
274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289
                  if (k .eq. kmax) then
                     vsd_v(i,j,k)=-dt*dry_v(i,j)*_HALF_*               &
                                   (tausy(i,j)+tausy(i,j+1))*rho_0i    &
                                  +auxo(k-1)*(vv(i,j,k)/hvo(i,j,k)     &
                                  -vv(i,j,k-1)/hvo(i,j,k-1))
                  end if
                  if ((k .gt. kvmin(i,j)) .and. (k .lt. kmax)) then
                     vsd_v(i,j,k)=-auxo(k)*(vv(i,j,k+1)/hvo(i,j,k+1)   &
                                  -vv(i,j,k)/hvo(i,j,k))               &
                                  +auxo(k-1)*(vv(i,j,k)/hvo(i,j,k)     &
                                  -vv(i,j,k-1)/hvo(i,j,k-1))
                  end if
                  if (k .eq. kvmin(i,j)) then
                     vsd_v(i,j,k)=-auxo(k)*(vv(i,j,k+1)/hvo(i,j,k+1)   &
                                  -vv(i,j,k)/hvo(i,j,k))
                  end if
Hans Burchard's avatar
Hans Burchard committed
290
#endif
291
#ifndef NO_BAROTROPIC
gotm's avatar
gotm committed
292
                  vv(i,j,k)=Res(k)+hvn(i,j,k)*Diff
293 294
#else
                  vv(i,j,k)=Res(k)
Hans Burchard's avatar
Hans Burchard committed
295 296 297
#endif
#ifdef _MOMENTUM_TERMS_
                  tdv_v(i,j,k)=(vv(i,j,k)-tdv_v(i,j,k))/dt
Hans Burchard's avatar
Hans Burchard committed
298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317
                  if (k .eq. kmax) then
                     vsd_v(i,j,k)=(vsd_v(i,j,k)                        &
                                  +auxn(k-1)*(vv(i,j,k)/hvn(i,j,k)     &
                                  -vv(i,j,k-1)/hvn(i,j,k-1)))/dt
                  end if
                  if ((k .gt. kvmin(i,j)) .and. (k .lt. kmax)) then
                     vsd_v(i,j,k)=(vsd_v(i,j,k)                        &
                                 -auxn(k)*(vv(i,j,k+1)/hvn(i,j,k+1)    &
                                  -vv(i,j,k)/hvn(i,j,k))               &
                                 +auxn(k-1)*(vv(i,j,k)/hvn(i,j,k)      &
                                  -vv(i,j,k-1)/hvn(i,j,k-1)))/dt
                  endif
                  if (k .eq. kvmin(i,j)) then
                     vsd_v(i,j,k)=(vsd_v(i,j,k)                        &
                                 -auxn(k)*(vv(i,j,k+1)/hvn(i,j,k+1)    &
                                  -vv(i,j,k)/hvn(i,j,k))               &
                                 +dt*rrv(i,j)*vv(i,j,k)                &
                                  /(_HALF_*(hvn(i,j,k)+hvo(i,j,k))))/dt
                  endif

318
#endif
gotm's avatar
gotm committed
319 320 321
               end do
            else ! (kmax .eq. kvmin(i,j))
               vv(i,j,kmax)=Vint(i,j)
kbk's avatar
kbk committed
322
            end if
gotm's avatar
gotm committed
323 324 325
         end if
      end do
   end do
bjb's avatar
bjb committed
326
!$OMP END DO
gotm's avatar
gotm committed
327

328
#ifdef SLICE_MODEL
bjb's avatar
bjb committed
329
!$OMP DO SCHEDULE(RUNTIME)
330
   do i=imin,imax
331 332 333 334 335
      do k=kvmin(i,2),kmax
         vv(i,1,k)=vv(i,2,k)
         vv(i,3,k)=vv(i,2,k)
      end do
   end do
bjb's avatar
bjb committed
336
!$OMP END DO
337 338
#endif

bjb's avatar
bjb committed
339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368
! Each thread must deallocate its own HEAP storage:
   deallocate(dif,stat=rc)
   if (rc /= 0) stop 'vv_momentum_3d: Error deallocating memory (dif)'

   deallocate(auxn,stat=rc)
   if (rc /= 0) stop 'vv_momentum_3d: Error deallocating memory (auxn)'

   deallocate(auxo,stat=rc)
   if (rc /= 0) stop 'vv_momentum_3d: Error deallocating memory (auxo)'

   deallocate(a1,stat=rc)
   if (rc /= 0) stop 'vv_momentum_3d: Error deallocating memory (a1)'

   deallocate(a2,stat=rc)
   if (rc /= 0) stop 'vv_momentum_3d: Error deallocating memory (a2)'

   deallocate(a3,stat=rc)
   if (rc /= 0) stop 'vv_momentum_3d: Error deallocating memory (a3)'

   deallocate(a4,stat=rc)
   if (rc /= 0) stop 'vv_momentum_3d: Error deallocating memory (a4)'

   deallocate(Res,stat=rc)
   if (rc /= 0) stop 'vv_momentum_3d: Error deallocating memory (Res)'

   deallocate(ex,stat=rc)
   if (rc /= 0) stop 'vv_momentum_3d: Error deallocating memory (ex)'

!$OMP END PARALLEL

gotm's avatar
gotm committed
369
!  Update the halo zones
bjb's avatar
bjb committed
370
   call tic(TIM_VVMOMENTUMH)
371
   call update_3d_halo(vv,vv,av,imin,jmin,imax,jmax,kmax,V_TAG)
gotm's avatar
gotm committed
372 373 374 375 376 377

   if (bdy3d) then
!      call do_bdy_3d(2,vv)
   end if

   call wait_halo(V_TAG)
bjb's avatar
bjb committed
378
   call toc(TIM_VVMOMENTUMH)
kbk's avatar
kbk committed
379
   call mirror_bdy_3d(vv,V_TAG)
gotm's avatar
gotm committed
380

381
   if (vel_check .ne. 0 .and. mod(n,abs(vel_check)) .eq. 0) then
382 383
      call check_3d_fields(imin,jmin,imax,jmax,kvmin,kmax,av, &
                           vv,min_vel,max_vel,status)
384 385 386 387 388 389
      if (status .gt. 0) then
         if (vel_check .gt. 0) then
            call getm_error("vv_momentum_3d()", &
                            "out-of-bound values encountered")
         end if
         if (vel_check .lt. 0) then
kbk's avatar
kbk committed
390
            LEVEL1 'vv_momentum_3d(): ',status, &
391 392 393 394 395
                   ' out-of-bound values encountered'
         end if
      end if
   end if

bjb's avatar
bjb committed
396
   call toc(TIM_VVMOMENTUM)
gotm's avatar
gotm committed
397 398 399 400 401 402 403 404 405 406 407
#ifdef DEBUG
   write(debug,*) 'Leaving vv_momentum_3d()'
   write(debug,*)
#endif
   return
   end subroutine vv_momentum_3d
!EOC

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