momentum.F90 11.8 KB
Newer Older
kbk's avatar
kbk committed
1
!$Id: momentum.F90,v 1.12 2006-03-01 15:54:07 kbk Exp $
gotm's avatar
gotm committed
2 3 4 5
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
6
! !IROUTINE: momentum - 2D-momentum for all interior points.
gotm's avatar
gotm committed
7 8 9 10 11 12
!
! !INTERFACE:
   subroutine momentum(n,tausx,tausy,airp)
!
! !DESCRIPTION:
!
hb's avatar
hb committed
13 14 15 16 17
! This small routine calls the $U$-equation and the $V$-equation in an
! alternating sequence (UVVUUVVUUVVU), in order to provide higher
! accuracy and energy conservation for the explicit numerical treatment 
! of the Coriolis term.
!
gotm's avatar
gotm committed
18 19 20 21 22
! !USES:
   use domain, only: imin,imax,jmin,jmax
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
kbk's avatar
kbk committed
23 24 25 26
   integer, intent(in)                 :: n
   REALTYPE, intent(in)                :: tausx(E2DFIELD)
   REALTYPE, intent(in)                :: tausy(E2DFIELD)
   REALTYPE, intent(in)                :: airp(E2DFIELD)
gotm's avatar
gotm committed
27 28 29 30 31
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
kbk's avatar
kbk committed
32 33 34
! !REVISION HISTORY:
!  Original author(s): Hans Burchard & Karsten Bolding
!
gotm's avatar
gotm committed
35
! !LOCAL VARIABLES:
kbk's avatar
kbk committed
36
  logical                    :: ufirst=.false.
gotm's avatar
gotm committed
37 38 39 40 41 42 43 44 45
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'Momentum() # ',Ncall
#endif

kbk's avatar
kbk committed
46
   if(ufirst) then
gotm's avatar
gotm committed
47 48
      call umomentum(tausx,airp)
      call vmomentum(tausy,airp)
kbk's avatar
kbk committed
49
      ufirst = .false.
gotm's avatar
gotm committed
50 51 52
   else
      call vmomentum(tausy,airp)
      call umomentum(tausx,airp)
kbk's avatar
kbk committed
53
      ufirst = .true.
gotm's avatar
gotm committed
54 55 56 57 58 59 60 61 62
   end if

   return
   end subroutine momentum
!EOC

!-----------------------------------------------------------------------
!BOP
!
63
! !IROUTINE: umomentum - 2D-momentum for all interior points.
gotm's avatar
gotm committed
64 65 66 67 68 69
!
! !INTERFACE:
   subroutine umomentum(tausx,airp)
!
! !DESCRIPTION:
!
hb's avatar
hb committed
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
! Here, the vertically integrated $U$-momentum equation (\ref{UMOM}) given
! on page \pageref{UMOM} including a
! number of slow terms is calculated. One slight modification is that
! for better stability of drying and flooding processes the slow friction
! term $S^x_F$ is now also multiplied with the parameter $\alpha$ defined
! in eq.\ (\ref{alpha}). 
!
! Furthermore, the horizontal pressure gradient $\partial^*_x\zeta$ is 
! modified in order to
! support drying and flooding, see figure \ref{figpressgrad} on page 
! \pageref{figpressgrad} and the explanations in section \ref{Section_dry}.
! $\partial^*_x\zeta$ is now also considering the atmospheric pressure
! gradient at sea surface height.
!
! For numerical stability reasons, the $U$-momentum equation is here 
! discretised in time such that the
! bed friction is treated explicitely:
!
!  \begin{equation}\label{Umom_discrete}
!  \displaystyle
!  U^{n+1}=\frac{U^n-\Delta t_m(gD\partial^*_x\zeta
!  +\alpha(-\frac{\tau_x^s}{\rho_0}-fV^n+U_{Ex}+S_A^x-S_D^x+S_B^x+S_F^x))}
!  {1+\Delta t_m\frac{R}{D^2}\sqrt{\left(U^n\right)^2+\left(V^n\right)^2}},
!  \end{equation}
!
!  with $U_{Ex}$ combining advection and diffusion of $U$, see routines
!  {\tt uv\_advect} (section \ref{sec-uv-advect} on page 
!  \pageref{sec-uv-advect}) and {\tt uv\_diffusion} 
!  (section \ref{sec-uv-diffusion} on page 
!  \pageref{sec-uv-diffusion}). The slow terms
!  are calculated in the routine {\tt slow\_terms} documented in section
!  \ref{sec-slow-terms} on page \pageref{sec-slow-terms}.
!  In (\ref{Umom_discrete}), $U^{n+1}$ denotes the transport on the
!  new and $U^n$ and $V^n$ the transports on the old time level.
!
!  The Coriolis term $fU$ for the subsequent $V$-momentum is also calculated
!  here, by directly interpolating the $U$-transports to the V-points
!  or by a method suggested by \cite{ESPELIDea00} which takes the
!  varying water depths into account.
!
!  Some provisions for proper behaviour of the $U$-transports when
!  GETM runs as slice model are made as well, see section
!  \ref{Section_GETM_Slice} on page \pageref{Section_GETM_Slice}.
!
gotm's avatar
gotm committed
114 115
! !USES:
   use parameters, only: g,rho_0
kbk's avatar
kbk committed
116 117
   use domain, only: imin,imax,jmin,jmax
   use domain, only: H,au,av,min_depth,dry_u,Cori,corv
gotm's avatar
gotm committed
118
#if defined(SPHERICAL) || defined(CURVILINEAR)
kbk's avatar
kbk committed
119
   use domain, only: dxu,arvd1,dxc,dyx
gotm's avatar
gotm committed
120 121
   use variables_2d, only: V
#else
kbk's avatar
kbk committed
122
   use domain, only: dx
gotm's avatar
gotm committed
123 124
#endif
   use m2d, only: dtm
kbk's avatar
kbk committed
125 126
   use variables_2d, only: D,z,UEx,U,DU,fV,SlUx,SlRu,ru,fU,DV
   use halo_zones, only : update_2d_halo,wait_halo,U_TAG
gotm's avatar
gotm committed
127 128 129
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
kbk's avatar
kbk committed
130
   REALTYPE, intent(in)                :: tausx(E2DFIELD),airp(E2DFIELD)
gotm's avatar
gotm committed
131 132 133 134 135 136
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !LOCAL VARIABLES:
kbk's avatar
kbk committed
137 138 139 140 141 142
   integer                   :: i,j
   REALTYPE                  :: zx(E2DFIELD)
   REALTYPE                  :: Slr(E2DFIELD),tausu(E2DFIELD)
   REALTYPE                  :: zp,zm,Uloc,Uold
   REALTYPE                  :: gamma=rho_0*g
   REALTYPE                  :: cord_curv=_ZERO_
gotm's avatar
gotm committed
143 144 145 146 147 148 149 150 151 152 153 154 155 156 157
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'umomentum() # ',Ncall
#endif

   do j=jmin,jmax
      do i=imin,imax
         if (au(i,j) .gt. 0) then
            zp=max(z(i+1,j),-H(i  ,j)+min(min_depth,D(i+1,j)))
            zm=max(z(i  ,j),-H(i+1,j)+min(min_depth,D(i  ,j)))
            zx(i,j)=(zp-zm+(airp(i+1,j)-airp(i,j))/gamma)/DXU
158
            tausu(i,j)=0.5*(tausx(i,j)+tausx(i+1,j))
gotm's avatar
gotm committed
159 160 161 162
         end if
      end do
   end do

kbk's avatar
kbk committed
163
   where (U .gt. 0)
gotm's avatar
gotm committed
164
      Slr=max(Slru, _ZERO_ )
kbk's avatar
kbk committed
165
   else where
gotm's avatar
gotm committed
166 167
      Slr=min(Slru, _ZERO_ )
   end where
168

kbk's avatar
kbk committed
169
   where ((au .eq. 1) .or. (au .eq. 2))
gotm's avatar
gotm committed
170 171 172
      U=(U-dtm*(g*DU*zx+dry_u*(-tausu/rho_0-fV+UEx+SlUx+Slr)))/(1+dtm*ru/DU)
   end where

173 174 175 176 177 178
#ifdef SLICE_MODEL
   do i=imin,imax
      U(i,3)=U(i,2)
   end do
#endif

gotm's avatar
gotm committed
179 180 181
!  now u is calculated
   call update_2d_halo(U,U,au,imin,jmin,imax,jmax,U_TAG)
   call wait_halo(U_TAG)
kbk's avatar
kbk committed
182
   call mirror_bdy_2d(U,U_TAG)
gotm's avatar
gotm committed
183 184 185 186

! Semi-implicit treatment of Coriolis force for V-momentum eq.
   do j=jmin,jmax
      do i=imin,imax
187
         if(av(i,j) .ge. 1) then
gotm's avatar
gotm committed
188 189
! Espelid et al. [2000], IJNME 49, 1521-1545
#ifdef NEW_CORI
kbk's avatar
kbk committed
190
            Uloc= &
kbk's avatar
kbk committed
191 192
             ( U(i,j  )/sqrt(DU(i,j  ))+ U(i-1,j  )/sqrt(DU(i-1,j  ))  &
             + U(i,j+1)/sqrt(DU(i,j+1))+ U(i-1,j+1)/sqrt(DU(i-1,j+1))) &
kbk's avatar
kbk committed
193
               *0.25*sqrt(DV(i,j))
gotm's avatar
gotm committed
194
#else
195
            Uloc=0.25*( U(i,j)+ U(i-1,j)+ U(i,j+1)+ U(i-1,j+1))
gotm's avatar
gotm committed
196 197
#endif
#if defined(SPHERICAL) || defined(CURVILINEAR)
kbk's avatar
kbk committed
198
            cord_curv=(V(i,j)*(DYX-DYXIM1)-Uloc*(DXCJP1-DXC)) &
kbk's avatar
kbk committed
199
                      /DV(i,j)*ARVD1
kbk's avatar
kbk committed
200
            fU(i,j)=(cord_curv+corv(i,j))*Uloc
gotm's avatar
gotm committed
201
#else
kbk's avatar
kbk committed
202
            fU(i,j)=corv(i,j)*Uloc
gotm's avatar
gotm committed
203
#endif
204
         else
kbk's avatar
kbk committed
205
            fU(i,j)= _ZERO_
206
         end if
gotm's avatar
gotm committed
207 208 209 210 211 212 213 214 215 216 217 218 219 220
      end do
   end do

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

!-----------------------------------------------------------------------
!BOP
!
221
! !IROUTINE: vmomentum - 2D-momentum for all interior points.
gotm's avatar
gotm committed
222 223 224 225 226 227
!
! !INTERFACE:
   subroutine vmomentum(tausy,airp)
!
! !DESCRIPTION:
!
hb's avatar
hb committed
228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271
! Here, the vertically integrated $V$-momentum equation (\ref{VMOM}) given
! on page \pageref{VMOM} including a
! number of slow terms is calculated. One slight modification is that
! for better stability of drying and flooding processes the slow friction
! term $S^y_F$ is now also multiplied with the parameter $\alpha$ defined
! in eq.\ (\ref{alpha}). 
!
! Furthermore, the horizontal pressure gradient $\partial^*_y\zeta$ is 
! modified in order to
! support drying and flooding, see figure \ref{figpressgrad} on page 
! \pageref{figpressgrad} and the explanations in section \ref{Section_dry}.
! $\partial^*_y\zeta$ is now also considering the atmospheric pressure
! gradient at sea surface height.
!
! For numerical stability reasons, the $V$-momentum equation is here 
! discretised in time such that the
! bed friction is treated explicitely:
!
!  \begin{equation}\label{Vmom_discrete}
!  \displaystyle
!  V^{n+1}=\frac{V^n-\Delta t_m(gD\partial^*_y\zeta
!  +\alpha(-\frac{\tau_y^s}{\rho_0}+fU^n+V_{Ex}+S_A^y-S_D^y+S_B^y+S_F^y))}
!  {1+\Delta t_m\frac{R}{D^2}\sqrt{\left(U^n\right)^2+\left(V^n\right)^2}},
!  \end{equation}
!
!  with $V_{Ex}$ combining advection and diffusion of $V$, see routines
!  {\tt uv\_advect} (section \ref{sec-uv-advect} on page
!  \pageref{sec-uv-advect}) and {\tt uv\_diffusion}
!  (section \ref{sec-uv-diffusion} on page
!  \pageref{sec-uv-diffusion}). The slow terms
!  are calculated in the routine {\tt slow\_terms} documented in section
!  \ref{sec-slow-terms} on page \pageref{sec-slow-terms}.
!  In (\ref{Vmom_discrete}), $V^{n+1}$ denotes the transport on the
!  new and $U^n$ and $V^n$ the transports on the old time level.
!
!  The Coriolis term $fV$ for the subsequent $U$-momentum is also calculated
!  here, by directly interpolating the $U$-transports to the U-points
!  or by a method suggested by \cite{ESPELIDea00} which takes the
!  varying water depths into account.
!
!  Some provisions for proper behaviour of the $V$-transports when
!  GETM runs as slice model are made as well, see section
!  \ref{Section_GETM_Slice} on page \pageref{Section_GETM_Slice}.
!
gotm's avatar
gotm committed
272 273
! !USES:
   use parameters, only: g,rho_0
kbk's avatar
kbk committed
274 275
   use domain, only: imin,imax,jmin,jmax
   use domain, only: H,au,av,min_depth,dry_v,Cori,coru
gotm's avatar
gotm committed
276
#if defined(SPHERICAL) || defined(CURVILINEAR)
kbk's avatar
kbk committed
277 278
   use domain, only: dyv,arud1,dxx,dyc
   use m2d, only: U
gotm's avatar
gotm committed
279
#else
kbk's avatar
kbk committed
280
   use domain, only: dy
gotm's avatar
gotm committed
281
#endif
kbk's avatar
kbk committed
282
   use m2d, only: dtm
kbk's avatar
kbk committed
283 284
   use variables_2d, only: D,z,VEx,V,DV,fU,SlVx,SlRv,rv,fV,DU
   use halo_zones, only : update_2d_halo,wait_halo,V_TAG
gotm's avatar
gotm committed
285 286 287
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
kbk's avatar
kbk committed
288
   REALTYPE, intent(in)                :: tausy(E2DFIELD),airp(E2DFIELD)
gotm's avatar
gotm committed
289 290 291 292 293 294
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !LOCAL VARIABLES:
kbk's avatar
kbk committed
295 296 297 298 299 300
   integer                   :: i,j
   REALTYPE                  :: zy(E2DFIELD)
   REALTYPE                  :: Slr(E2DFIELD),tausv(E2DFIELD)
   REALTYPE                  :: zp,zm,Vloc
   REALTYPE                  :: gamma=rho_0*g
   REALTYPE                  :: cord_curv=_ZERO_
gotm's avatar
gotm committed
301 302 303 304 305 306 307 308 309 310 311 312 313 314 315
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'vmomentum() # ',Ncall
#endif

   do j=jmin,jmax
      do i=imin,imax
         if (av(i,j) .gt. 0) then
            zp=max(z(i,j+1),-H(i,j  )+min(min_depth,D(i,j+1)))
            zm=max(z(i,j  ),-H(i,j+1)+min(min_depth,D(i,j  )))
            zy(i,j)=(zp-zm+(airp(i,j+1)-airp(i,j))/gamma)/DYV
316
            tausv(i,j)=0.5*(tausy(i,j)+tausy(i,j+1))
gotm's avatar
gotm committed
317 318 319 320
         end if
      end do
   end do

kbk's avatar
kbk committed
321
   where (V .gt. 0)
gotm's avatar
gotm committed
322
      Slr=max(Slrv, _ZERO_ )
kbk's avatar
kbk committed
323
   else where
gotm's avatar
gotm committed
324 325
      Slr=min(Slrv, _ZERO_ )
   end where
326

kbk's avatar
kbk committed
327
   where ((av .eq. 1) .or. (av .eq. 2))
gotm's avatar
gotm committed
328 329 330
      V=(V-dtm*(g*DV*zy+dry_v*(-tausv/rho_0+fU+VEx+SlVx+Slr)))/(1+dtm*rv/DV)
   end where

331 332 333 334 335 336 337
#ifdef SLICE_MODEL
   do i=imin,imax
      V(i,1)=V(i,2)
      V(i,3)=V(i,2)
   end do
#endif

gotm's avatar
gotm committed
338 339 340
!  now v is calculated
   call update_2d_halo(V,V,av,imin,jmin,imax,jmax,V_TAG)
   call wait_halo(V_TAG)
kbk's avatar
kbk committed
341
   call mirror_bdy_2d(V,V_TAG)
gotm's avatar
gotm committed
342 343 344 345

!  Semi-implicit treatment of Coriolis force for U-momentum eq.
   do j=jmin,jmax
      do i=imin,imax
346
         if(au(i,j) .ge. 1) then
gotm's avatar
gotm committed
347 348
! Espelid et al. [2000], IJNME 49, 1521-1545
#ifdef NEW_CORI
kbk's avatar
kbk committed
349
            Vloc = &
kbk's avatar
kbk committed
350 351 352
            ( V(i,j  )/sqrt(DV(i,j  ))+ V(i+1,j  )/sqrt(DV(i+1,j  )) + &
              V(i,j-1)/sqrt(DV(i,j-1))+ V(i+1,j-1)/sqrt(DV(i+1,j-1)))  &
              *0.25*sqrt(DU(i,j))
gotm's avatar
gotm committed
353
#else
kbk's avatar
kbk committed
354
            Vloc = 0.25*( V(i,j)+ V(i+1,j)+ V(i,j-1)+ V(i+1,j-1))
gotm's avatar
gotm committed
355 356
#endif
#if defined(SPHERICAL) || defined(CURVILINEAR)
kbk's avatar
kbk committed
357
            cord_curv=(Vloc*(DYCIP1-DYC)-U(i,j)*(DXX-DXXJM1)) &
kbk's avatar
kbk committed
358
                       /DU(i,j)*ARUD1
kbk's avatar
kbk committed
359
            fV(i,j)=(cord_curv+coru(i,j))*Vloc
gotm's avatar
gotm committed
360
#else
kbk's avatar
kbk committed
361
            fV(i,j)=coru(i,j)*Vloc
gotm's avatar
gotm committed
362
#endif
363 364 365
         else
            fV(i,j) = _ZERO_
         end if
gotm's avatar
gotm committed
366 367 368 369 370 371 372 373 374 375 376 377 378 379
      end do
   end do

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

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