New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
wzvmod.F90 in trunk/NEMO/OPA_SRC/DYN – NEMO

source: trunk/NEMO/OPA_SRC/DYN/wzvmod.F90 @ 708

Last change on this file since 708 was 708, checked in by smasson, 17 years ago

continue changeset:704, see ticket:5
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 8.2 KB
Line 
1MODULE wzvmod
2   !!==============================================================================
3   !!                       ***  MODULE  wzvmod  ***
4   !! Ocean diagnostic variable : vertical velocity
5   !!==============================================================================
6   !! History :   5.0  !  90-10  (C. Levy, G. Madec)  Original code
7   !!             7.0  !  96-01  (G. Madec)  Statement function for e3
8   !!             8.5  !  02-07  (G. Madec)  Free form, F90
9   !!----------------------------------------------------------------------
10   !!   wzv        : Compute the vertical velocity
11   !!----------------------------------------------------------------------
12   !! * Modules used
13   USE oce             ! ocean dynamics and tracers variables
14   USE dom_oce         ! ocean space and time domain variables
15   USE sbc_oce         ! surface boundary condition: ocean
16   USE domvvl          ! Variable volume
17   USE in_out_manager  ! I/O manager
18   USE prtctl          ! Print control
19   USE phycst
20   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
21
22   IMPLICIT NONE
23   PRIVATE
24
25   !! * Routine accessibility
26   PUBLIC wzv          ! routine called by step.F90 and inidtr.F90
27
28   !! * Substitutions
29#  include "domzgr_substitute.h90"
30   !!----------------------------------------------------------------------
31   !!  OPA 9.0 , LOCEAN-IPSL (2005)
32   !! $Id$
33   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
34   !!----------------------------------------------------------------------
35
36CONTAINS
37
38#if defined key_mpp_omp
39   !!----------------------------------------------------------------------
40   !!   'key_mpp_omp'                                   j-k-i loop (j-slab)
41   !!----------------------------------------------------------------------
42
43   SUBROUTINE wzv( kt )
44      !!----------------------------------------------------------------------
45      !!                    ***  ROUTINE wzv  ***
46      !!                     
47      !! ** Purpose :   Compute the now vertical velocity after the array swap
48      !!
49      !! ** Method  :   Using the incompressibility hypothesis, the vertical
50      !!     velocity is computed by integrating the horizontal divergence
51      !!     from the bottom to the surface.
52      !!     The boundary conditions are w=0 at the bottom (no flux) and,
53      !!     in rigid-lid case, w=0 at the sea surface.
54      !!
55      !! ** action  :    wn array : the now vertical velocity
56      !!----------------------------------------------------------------------
57      !! * Arguments
58      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
59
60      !! * Local declarations
61      INTEGER ::   jj, jk      ! dummy loop indices
62      !!----------------------------------------------------------------------
63
64      IF( kt == nit000 ) THEN
65         IF(lwp) WRITE(numout,*)
66         IF(lwp) WRITE(numout,*) 'wzv     : vertical velocity from continuity eq.'
67         IF(lwp) WRITE(numout,*) '~~~~~~~                     j-k-i loops'
68
69         ! bottom boundary condition: w=0 (set once for all)
70         wn(:,:,jpk) = 0.e0
71      ENDIF
72
73      !                                                ! ===============
74      DO jj = 1, jpj                                   !  Vertical slab
75         !                                             ! ===============
76         ! Computation from the bottom
77         DO jk = jpkm1, 1, -1
78            wn(:,jj,jk) = wn(:,jj,jk+1) - fse3t(:,jj,jk) * hdivn(:,jj,jk)
79         END DO
80         !                                             ! ===============
81      END DO                                           !   End of slab
82      !                                                ! ===============
83
84      IF(ln_ctl)   CALL prt_ctl(tab3d_1=wn, clinfo1=' w**2 -   : ', mask1=wn)
85
86   END SUBROUTINE wzv
87
88#else
89   !!----------------------------------------------------------------------
90   !!   Default option                                           k-j-i loop
91   !!----------------------------------------------------------------------
92
93   SUBROUTINE wzv( kt )
94      !!----------------------------------------------------------------------
95      !!                    ***  ROUTINE wzv  ***
96      !!
97      !! ** Purpose :   Compute the now vertical velocity after the array swap
98      !!
99      !! ** Method  :   Using the incompressibility hypothesis, the vertical
100      !!      velocity is computed by integrating the horizontal divergence
101      !!      from the bottom to the surface.
102      !!        The boundary conditions are w=0 at the bottom (no flux) and,
103      !!      in regid-lid case, w=0 at the sea surface.
104      !!
105      !! ** action  :   wn array : the now vertical velocity
106      !!----------------------------------------------------------------------
107      !! * Arguments
108      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
109
110      !! * Local declarations
111      INTEGER  ::           jk           ! dummy loop indices
112      !! Variable volume
113      INTEGER  ::   ji, jj               ! dummy loop indices
114      REAL(wp) ::   z2dt, zraur          ! temporary scalar
115      REAL(wp), DIMENSION (jpi,jpj) ::   zssha, zun, zvn, zhdiv
116      !!----------------------------------------------------------------------
117
118      IF( kt == nit000 ) THEN
119         IF(lwp) WRITE(numout,*)
120         IF(lwp) WRITE(numout,*) 'wzv     : vertical velocity from continuity eq.'
121         IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
122
123         ! bottom boundary condition: w=0 (set once for all)
124         wn(:,:,jpk) = 0.e0
125      ENDIF
126
127      IF( lk_vvl ) THEN                ! Variable volume
128         !
129         z2dt = 2. * rdt                                       ! time step: leap-frog
130         IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt       ! time step: Euler if restart from rest
131         zraur  = 1. / rauw
132
133         ! Vertically integrated quantities
134         ! --------------------------------
135         zun(:,:) = 0.e0
136         zvn(:,:) = 0.e0
137         !
138         DO jk = 1, jpkm1             ! Vertically integrated transports (now)
139            zun(:,:) = zun(:,:) + fse3u(:,:,jk) * un(:,:,jk)
140            zvn(:,:) = zvn(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
141         END DO
142
143         ! Horizontal divergence of barotropic transports
144         !--------------------------------------------------
145         DO jj = 2, jpjm1
146            DO ji = 2, jpim1   ! vector opt.
147               zhdiv(ji,jj) = (  e2u(ji  ,jj  ) * zun(ji  ,jj  )     &
148                  &            - e2u(ji-1,jj  ) * zun(ji-1,jj  )     &
149                  &            + e1v(ji  ,jj  ) * zvn(ji  ,jj  )     &
150                  &            - e1v(ji  ,jj-1) * zvn(ji  ,jj-1) )   &
151                  &           / ( e1t(ji,jj) * e2t(ji,jj) )
152            END DO
153         END DO
154
155#if defined key_obc && ( key_dynspg_exp || key_dynspg_ts )
156         ! open boundaries (div must be zero behind the open boundary)
157         !  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
158         IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1)   = 0.e0    ! east
159         IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1)   = 0.e0    ! west
160         IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0    ! north
161         IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1)   = 0.e0    ! south
162#endif
163
164         CALL lbc_lnk( zhdiv, 'T', 1. )
165
166         ! Sea surface elevation time stepping
167         ! -----------------------------------
168         zssha(:,:) = sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) ) * tmask(:,:,1)
169
170         ! Vertical velocity computed from bottom
171         ! --------------------------------------
172         DO jk = jpkm1, 1, -1
173            wn(:,:,jk) = wn(:,:,jk+1) - fse3t(:,:,jk) * hdivn(:,:,jk) &
174              &        - ( zssha(:,:) - sshb(:,:) ) * fsve3t(:,:,jk) * mut(:,:,jk) / z2dt
175         END DO
176
177      ELSE                             ! Fixed volume
178
179         ! Vertical velocity computed from bottom
180         ! --------------------------------------
181         DO jk = jpkm1, 1, -1
182            wn(:,:,jk) = wn(:,:,jk+1) - fse3t(:,:,jk) * hdivn(:,:,jk)
183         END DO
184     
185      ENDIF
186
187      IF(ln_ctl)   CALL prt_ctl(tab3d_1=wn, clinfo1=' w**2 -   : ', mask1=wn)
188
189   END SUBROUTINE wzv
190#endif
191
192   !!======================================================================
193END MODULE wzvmod
Note: See TracBrowser for help on using the repository browser.