uu_momentum_3d.F90 14.4 KB
Newer Older
gotm's avatar
gotm committed
1 2 3 4
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
hb's avatar
hb committed
5
! !ROUTINE: uu_momentum_3d - $x$-momentum eq.\ \label{sec-uu-momentum-3d}
gotm's avatar
gotm committed
6 7
!
! !INTERFACE:
8
   subroutine uu_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 21 22 23 24 25 26 27
! Here, the budget equation for layer-averaged momentum in eastern direction,
! $p_k$,
! is calculated. The physical equation is given as equation (\ref{uEq}),
! the layer-integrated equation as (\ref{uEqvi}), and after curvilinear
! transformation as (\ref{uEqviCurvi}).
! In this routine, first the Coriolis rotation term, $fq_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).
!
! 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 U-point,
! and the barotropic pressure gradient is calculated (the latter
! includes the pressure gradient correction for drying points, see
! section \ref{Section_dry}).
28 29
! Afterwards, the matrix is set up for each water column, and it is solved
! by means of a tri-diagonal matrix solver.
hb's avatar
hb committed
30
!
31
! In case that the compiler option {\tt STRUCTURE\_FRICTION} is switched on,
32 33 34 35 36 37
! the frictional effect of structures in the water column is calculated
! by adding the quadratic frictional term $C u \sqrt{u^2+v^2}$ (with a minus sign on
! the right hand side) numerically implicitly to the $u$-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}.
!
hb's avatar
hb committed
38 39 40
! Finally, the new velocity profile is shifted such that its vertical
! integral is identical to the time integral of the vertically integrated
! transport.
41
!
hb's avatar
hb committed
42
! When GETM is run as a slice model (compiler option {\tt SLICE\_MODEL}
43
! is activated), the result for $j=2$ is copied to $j=3$.
hb's avatar
hb committed
44
!
gotm's avatar
gotm committed
45
! !USES:
46
   use exceptions
gotm's avatar
gotm committed
47
   use parameters, only: g,avmmol,rho_0
Knut's avatar
Knut committed
48
   use domain, only: imin,imax,jmin,jmax,kmax,H,min_depth
49
   use domain, only: dry_u,coru,au,av,az,ax
gotm's avatar
gotm committed
50
#if defined CURVILINEAR || defined SPHERICAL
51
   use domain, only: dxu,arud1,dxx,dyc,dyx,dxc
gotm's avatar
gotm committed
52
#else
53
   use domain, only: dx,dy
gotm's avatar
gotm committed
54
#endif
55 56
   use variables_2d, only: Uint,D
   use bdy_3d, only: do_bdy_3d
gotm's avatar
gotm committed
57
   use variables_3d, only: dt,cnpar,kumin,uu,vv,huo,hun,hvo,uuEx,ww,hvn
Knut's avatar
Knut committed
58
   use variables_3d, only: num,nuh,sseo,Dun,rru
Hans Burchard's avatar
Hans Burchard committed
59
#ifdef _MOMENTUM_TERMS_
60
   use variables_3d, only: tdv_u,cor_u,ipg_u,epg_u,vsd_u,hsd_u
Hans Burchard's avatar
Hans Burchard committed
61
#endif
62 63 64
#ifdef STRUCTURE_FRICTION
   use variables_3d, only: sf
#endif
65 66
#ifndef NO_BAROCLINIC
   use variables_3d, only: idpdx
gotm's avatar
gotm committed
67
#endif
68 69
   use halo_zones, only: update_3d_halo,wait_halo,U_TAG
   use meteo, only: tausx,airp
kbk's avatar
kbk committed
70
   use m3d, only: ip_fac
71
   use m3d, only: vel_check,min_vel,max_vel
bjb's avatar
bjb committed
72
   use getm_timers, only: tic, toc, TIM_UUMOMENTUM, TIM_UUMOMENTUMH
bjb's avatar
bjb committed
73
!$ use omp_lib
gotm's avatar
gotm committed
74 75 76
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
77
   integer, intent(in)                 :: n
kbk's avatar
kbk committed
78
   logical, intent(in)                 :: bdy3d
gotm's avatar
gotm committed
79
!
kbk's avatar
kbk committed
80 81 82
! !REVISION HISTORY:
!  Original author(s): Hans Burchard & Karsten Bolding
!
gotm's avatar
gotm committed
83
! !LOCAL VARIABLES:
kbk's avatar
kbk committed
84
   integer                   :: i,j,k,rc
bjb's avatar
bjb committed
85 86 87 88 89
   REALTYPE, POINTER         :: dif(:)
   REALTYPE, POINTER         :: auxn(:),auxo(:)
   REALTYPE, POINTER         :: a1(:),a2(:)
   REALTYPE, POINTER         :: a3(:),a4(:)
   REALTYPE, POINTER         :: Res(:),ex(:)
kbk's avatar
kbk committed
90 91 92
   REALTYPE                  :: zp,zm,zx,ResInt,Diff,Vloc
   REALTYPE                  :: gamma=g*rho_0
   REALTYPE                  :: cord_curv=_ZERO_
bjb's avatar
bjb committed
93
   REALTYPE                  :: gammai,rho_0i
94
   integer                   :: status
gotm's avatar
gotm committed
95 96 97 98 99 100 101 102
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'uu_momentum_3d() # ',Ncall
#endif
bjb's avatar
bjb committed
103
   call tic(TIM_UUMOMENTUM)
104

bjb's avatar
bjb committed
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
   gammai=_ONE_/gamma
   rho_0i=_ONE_/rho_0
!$OMP PARALLEL DEFAULT(SHARED)                                         &
!$OMP    PRIVATE(i,j,k,rc,zp,zm,zx,ResInt,Diff,Vloc,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 'uu_momentum_3d: Error allocating memory (dif)'

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

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

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

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

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

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

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

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

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

!$OMP DO SCHEDULE(RUNTIME)
143 144
   do j=jmin,jmax
      do i=imin,imax
gotm's avatar
gotm committed
145 146

         if (au(i,j) .eq. 1 .or. au(i,j) .eq. 2) then
147
            if (kmax .gt. kumin(i,j)) then
gotm's avatar
gotm committed
148 149 150
               do k=kumin(i,j),kmax ! explicit terms
! Espelid et al. [2000], IJNME 49, 1521-1545
#ifdef NEW_CORI
151 152 153 154
                  Vloc=(vv(i  ,j  ,k)/sqrt(hvo(i  ,j  ,k))   &
                       +vv(i+1,j  ,k)/sqrt(hvo(i+1,j  ,k))   &
                       +vv(i  ,j-1,k)/sqrt(hvo(i  ,j-1,k))   &
                       +vv(i+1,j-1,k)/sqrt(hvo(i+1,j-1,k)))  &
bjb's avatar
bjb committed
155
                       *_QUART_*sqrt(huo(i,j,k))
gotm's avatar
gotm committed
156
#else
bjb's avatar
bjb committed
157
                  Vloc=_QUART_*(vv(i,j,k)+vv(i+1,j,k)+vv(i,j-1,k)+vv(i+1,j-1,k))
gotm's avatar
gotm committed
158
#endif
159 160
#if defined(SPHERICAL) || defined(CURVILINEAR)
                  cord_curv=(Vloc*(DYCIP1-DYC)-uu(i,j,k)*(DXX-DXXJM1))   &
gotm's avatar
gotm committed
161
                        /huo(i,j,k)*ARUD1
162
                  ex(k)=(cord_curv+coru(i,j))*Vloc
gotm's avatar
gotm committed
163
#else
164 165
                  ex(k)=coru(i,j)*Vloc
#endif
Hans Burchard's avatar
Hans Burchard committed
166 167 168
#ifdef _MOMENTUM_TERMS_
                  cor_u(i,j,k)=-dry_u(i,j)*ex(k)
#endif
169 170 171
#ifdef NO_BAROCLINIC
                  ex(k)=dry_u(i,j)*(ex(k)-uuEx(i,j,k))
#else
kbk's avatar
kbk committed
172
                  ex(k)=dry_u(i,j)*(ex(k)-uuEx(i,j,k)+ip_fac*idpdx(i,j,k))
Hans Burchard's avatar
Hans Burchard committed
173 174 175
#ifdef _MOMENTUM_TERMS_
                  ipg_u(i,j,k)=-dry_u(i,j)*ip_fac*idpdx(i,j,k)
#endif
gotm's avatar
gotm committed
176 177 178
#endif
               end do
               ex(kmax)=ex(kmax)                                         &
bjb's avatar
bjb committed
179
                       +dry_u(i,j)*_HALF_*(tausx(i,j)+tausx(i+1,j))*rho_0i
gotm's avatar
gotm committed
180 181 182 183 184 185 186
!     Eddy viscosity
               do k=kumin(i,j),kmax-1
                  dif(k)=0.5*(num(i,j,k)+num(i+1,j,k)) + avmmol
               end do

!     Auxilury terms, old and new time level,
               do k=kumin(i,j),kmax-1
bjb's avatar
bjb committed
187 188
                  auxo(k)=_TWO_*(_ONE_-cnpar)*dt*dif(k)/(huo(i,j,k+1)+huo(i,j,k))
                  auxn(k)=_TWO_*       cnpar *dt*dif(k)/(hun(i,j,k+1)+hun(i,j,k))
gotm's avatar
gotm committed
189 190 191
               end do

!     Barotropic pressure gradient
192
#ifndef NO_BAROTROPIC
gotm's avatar
gotm committed
193 194
               zp=max(sseo(i+1,j),-H(i  ,j)+min(min_depth,D(i+1,j)))
               zm=max(sseo(i  ,j),-H(i+1,j)+min(min_depth,D(i  ,j)))
bjb's avatar
bjb committed
195
               zx=(zp-zm+(airp(i+1,j)-airp(i,j))*gammai)/DXU
196
#else
bjb's avatar
bjb committed
197
               zx=_ZERO_
198
#endif
gotm's avatar
gotm committed
199 200 201 202 203

!     Matrix elements for surface layer
               k=kmax
               a1(k)=-auxn(k-1)/hun(i,j,k-1)
               a2(k)=1.+auxn(k-1)/hun(i,j,k)
kbk's avatar
kbk committed
204 205 206
               a4(k)=uu(i,j,k  )*(1-auxo(k-1)/huo(i,j,k))              &
                    +uu(i,j,k-1)*auxo(k-1)/huo(i,j,k-1)                &
                    +dt*ex(k)                                          &
bjb's avatar
bjb committed
207
                    -dt*_HALF_*(huo(i,j,k)+hun(i,j,k))*g*zx
gotm's avatar
gotm committed
208 209 210 211 212

!     Matrix elements for inner layers
               do k=kumin(i,j)+1,kmax-1
                  a3(k)=-auxn(k  )/hun(i,j,k+1)
                  a1(k)=-auxn(k-1)/hun(i,j,k-1)
bjb's avatar
bjb committed
213
                  a2(k)=_ONE_+(auxn(k)+auxn(k-1))/hun(i,j,k)
kbk's avatar
kbk committed
214 215 216
                  a4(k)=uu(i,j,k+1)*auxo(k)/huo(i,j,k+1)               &
                       +uu(i,j,k  )*(1-(auxo(k)+auxo(k-1))/huo(i,j,k)) &
                       +uu(i,j,k-1)*auxo(k-1)/huo(i,j,k-1)             &
217
                       +dt*ex(k)                                       &
bjb's avatar
bjb committed
218
                       -dt*_HALF_*(huo(i,j,k)+hun(i,j,k))*g*zx
gotm's avatar
gotm committed
219 220 221 222 223
               end do

!     Matrix elements for bottom layer
               k=kumin(i,j)
               a3(k)=-auxn(k  )/hun(i,j,k+1)
bjb's avatar
bjb committed
224 225
               a2(k)=_ONE_+auxn(k)/hun(i,j,k)                          &
                     +dt*rru(i,j)/(_HALF_*(hun(i,j,k)+huo(i,j,k)))
kbk's avatar
kbk committed
226 227 228
               a4(k)=uu(i,j,k+1)*auxo(k)/huo(i,j,k+1)                  &
                    +uu(i,j,k  )*(1-auxo(k)/huo(i,j,k))                &
                    +dt*ex(k)                                          &
bjb's avatar
bjb committed
229
                    -dt*_HALF_*(huo(i,j,k)+hun(i,j,k))*g*zx
gotm's avatar
gotm committed
230

231 232
#ifdef STRUCTURE_FRICTION
               do k=kumin(i,j),kmax
bjb's avatar
bjb committed
233
                  a2(k)=a2(k)+dt*_HALF_*(sf(i,j,k)+sf(i+1,j,k))
234 235 236
               end do
#endif

gotm's avatar
gotm committed
237 238 239 240 241 242 243 244 245
               call getm_tridiagonal(kmax,kumin(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, Uint.

               ResInt= _ZERO_
               do k=kumin(i,j),kmax
                  ResInt=ResInt+Res(k)
               end do
Knut's avatar
Knut committed
246
               Diff=(Uint(i,j)-ResInt)/Dun(i,j)
gotm's avatar
gotm committed
247 248

               do k=kumin(i,j),kmax
Hans Burchard's avatar
Hans Burchard committed
249 250 251 252
#ifdef _MOMENTUM_TERMS_
                  tdv_u(i,j,k)=uu(i,j,k)
                  epg_u(i,j,k)=_HALF_*(huo(i,j,k)+hun(i,j,k))*g*zx     &
                              -hun(i,j,k)*Diff/dt
253
                  hsd_u(i,j,k)=hsd_u(i,j,k)*dry_u(i,j)
Hans Burchard's avatar
Hans Burchard committed
254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269
                  if (k .eq. kmax) then
                     vsd_u(i,j,k)=-dt*dry_u(i,j)*_HALF_*               &
                                   (tausx(i,j)+tausx(i+1,j))*rho_0i    &
                                  +auxo(k-1)*(uu(i,j,k)/huo(i,j,k)     &
                                  -uu(i,j,k-1)/huo(i,j,k-1))
                  end if
                  if ((k .gt. kumin(i,j)) .and. (k .lt. kmax)) then
                     vsd_u(i,j,k)=-auxo(k)*(uu(i,j,k+1)/huo(i,j,k+1)   &
                                  -uu(i,j,k)/huo(i,j,k))               &
                                  +auxo(k-1)*(uu(i,j,k)/huo(i,j,k)     &
                                  -uu(i,j,k-1)/huo(i,j,k-1))
                  end if
                  if (k .eq. kumin(i,j)) then
                     vsd_u(i,j,k)=-auxo(k)*(uu(i,j,k+1)/huo(i,j,k+1)   &
                                  -uu(i,j,k)/huo(i,j,k))
                  end if
Hans Burchard's avatar
Hans Burchard committed
270
#endif
271
#ifndef NO_BAROTROPIC
gotm's avatar
gotm committed
272
                  uu(i,j,k)=Res(k) +hun(i,j,k)*Diff
273 274
#else
                  uu(i,j,k)=Res(k)
Hans Burchard's avatar
Hans Burchard committed
275 276 277
#endif
#ifdef _MOMENTUM_TERMS_
                  tdv_u(i,j,k)=(uu(i,j,k)-tdv_u(i,j,k))/dt
Hans Burchard's avatar
Hans Burchard committed
278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296
                  if (k .eq. kmax) then
                     vsd_u(i,j,k)=(vsd_u(i,j,k)                        &
                                  +auxn(k-1)*(uu(i,j,k)/hun(i,j,k)     &
                                  -uu(i,j,k-1)/hun(i,j,k-1)))/dt
                  end if
                  if ((k .gt. kumin(i,j)) .and. (k .lt. kmax)) then
                     vsd_u(i,j,k)=(vsd_u(i,j,k)                        &
                                 -auxn(k)*(uu(i,j,k+1)/hun(i,j,k+1)    &
                                  -uu(i,j,k)/hun(i,j,k))               &
                                 +auxn(k-1)*(uu(i,j,k)/hun(i,j,k)      &
                                  -uu(i,j,k-1)/hun(i,j,k-1)))/dt
                  endif
                  if (k .eq. kumin(i,j)) then
                     vsd_u(i,j,k)=(vsd_u(i,j,k)                        &
                                 -auxn(k)*(uu(i,j,k+1)/hun(i,j,k+1)    &
                                  -uu(i,j,k)/hun(i,j,k))               &
                                 +dt*rru(i,j)*uu(i,j,k)                &
                                  /(_HALF_*(hun(i,j,k)+huo(i,j,k))))/dt
                  endif
297
#endif
gotm's avatar
gotm committed
298 299 300 301 302 303
               end do
            else  ! if (kmax .eq. kumin(i,j))
                  uu(i,j,kmax)=Uint(i,j)
            end if
         end if
      end do
Hans Burchard's avatar
Hans Burchard committed
304
   end do
bjb's avatar
bjb committed
305
!$OMP END DO
306 307

#ifdef SLICE_MODEL
bjb's avatar
bjb committed
308
!$OMP DO SCHEDULE(RUNTIME)
309
   do i=imin,imax
310 311 312
      do k=kumin(i,2),kmax
         uu(i,3,k)=uu(i,2,k)
      end do
gotm's avatar
gotm committed
313
   end do
bjb's avatar
bjb committed
314
!$OMP END DO NOWAIT
315 316
#endif

bjb's avatar
bjb committed
317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345
! Each thread must deallocate its own HEAP storage:
   deallocate(dif,stat=rc)
   if (rc /= 0) stop 'uu_momentum_3d: Error deallocating memory (dif)'

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

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

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

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

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

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

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

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

!$OMP END PARALLEL
gotm's avatar
gotm committed
346 347

!  Update the halo zones
bjb's avatar
bjb committed
348
   call tic(TIM_UUMOMENTUMH)
349
   call update_3d_halo(uu,uu,au,imin,jmin,imax,jmax,kmax,U_TAG)
gotm's avatar
gotm committed
350 351 352 353 354 355

   if (bdy3d) then
!      call do_bdy_3d(1,uu)
   end if

   call wait_halo(U_TAG)
bjb's avatar
bjb committed
356
   call toc(TIM_UUMOMENTUMH)
kbk's avatar
kbk committed
357
   call mirror_bdy_3d(uu,U_TAG)
gotm's avatar
gotm committed
358

359
   if (vel_check .ne. 0 .and. mod(n,abs(vel_check)) .eq. 0) then
360 361
      call check_3d_fields(imin,jmin,imax,jmax,kumin,kmax,au, &
                           uu,min_vel,max_vel,status)
362 363 364 365 366 367
      if (status .gt. 0) then
         if (vel_check .gt. 0) then
            call getm_error("uu_momentum_3d()", &
                            "out-of-bound values encountered")
         end if
         if (vel_check .lt. 0) then
kbk's avatar
kbk committed
368
            LEVEL1 'uu_momentum_3d(): ',status, &
369 370 371 372 373
                   ' out-of-bound values encountered'
         end if
      end if
   end if

bjb's avatar
bjb committed
374
   call toc(TIM_UUMOMENTUM)
gotm's avatar
gotm committed
375 376 377 378 379 380 381 382 383 384 385
#ifdef DEBUG
   write(debug,*) 'Leaving uu_momentum_3d()'
   write(debug,*)
#endif
   return
   end subroutine uu_momentum_3d
!EOC

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