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 @ 592

Last change on this file since 592 was 592, checked in by opalod, 17 years ago

nemo_v2_update_001 : CT : - add non linear free surface (variable volume) with new cpp key key_vvl

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