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

Last change on this file since 1241 was 1241, checked in by rblod, 15 years ago

Fix a stupid bug for time splitting and ensure restartability for dynspg_ts in addition, see tickets #280 and #292
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 6.3 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   !!              "   !  07-07  (D. Storkey) Zero zhdiv at open boundary (BDY)
10   !!----------------------------------------------------------------------
11   !!   wzv        : Compute the vertical velocity
12   !!----------------------------------------------------------------------
13   !! * Modules used
14   USE oce             ! ocean dynamics and tracers variables
15   USE dom_oce         ! ocean space and time domain variables
16   USE sbc_oce         ! surface boundary condition: ocean
17   USE domvvl          ! Variable volume
18   USE in_out_manager  ! I/O manager
19   USE prtctl          ! Print control
20   USE phycst
21   USE bdy_oce         ! unstructured open boundaries
22   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
23   USE obc_par         ! open boundary cond. parameter
24   USE obc_oce
25
26   IMPLICIT NONE
27   PRIVATE
28
29   !! * Routine accessibility
30   PUBLIC wzv          ! routine called by step.F90 and inidtr.F90
31
32   !! * Substitutions
33#  include "domzgr_substitute.h90"
34   !!----------------------------------------------------------------------
35   !!  OPA 9.0 , LOCEAN-IPSL (2005)
36   !! $Id$
37   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
38   !!----------------------------------------------------------------------
39
40CONTAINS
41
42   SUBROUTINE wzv( kt )
43      !!----------------------------------------------------------------------
44      !!                    ***  ROUTINE wzv  ***
45      !!
46      !! ** Purpose :   Compute the now vertical velocity after the array swap
47      !!
48      !! ** Method  :   Using the incompressibility hypothesis, the vertical
49      !!      velocity is computed by integrating the horizontal divergence
50      !!      from the bottom to the surface.
51      !!        The boundary conditions are w=0 at the bottom (no flux) and,
52      !!      in regid-lid case, w=0 at the sea surface.
53      !!
54      !! ** action  :   wn array : the now vertical velocity
55      !!----------------------------------------------------------------------
56      !! * Arguments
57      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
58
59      !! * Local declarations
60      INTEGER  ::           jk           ! dummy loop indices
61      !! Variable volume
62      INTEGER  ::   ji, jj               ! dummy loop indices
63      REAL(wp) ::   z2dt, zraur          ! temporary scalar
64      REAL(wp), DIMENSION (jpi,jpj) ::   zssha, zun, zvn, zhdiv
65#if defined key_bdy
66      INTEGER  ::     jgrd, jb           ! temporary scalars
67#endif
68      !!----------------------------------------------------------------------
69
70      IF( kt == nit000 ) THEN
71         IF(lwp) WRITE(numout,*)
72         IF(lwp) WRITE(numout,*) 'wzv     : vertical velocity from continuity eq.'
73         IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
74
75         ! bottom boundary condition: w=0 (set once for all)
76         wn(:,:,jpk) = 0.e0
77      ENDIF
78
79      IF( lk_vvl ) THEN                ! Variable volume
80         !
81         z2dt = 2. * rdt                                       ! time step: leap-frog
82         IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt       ! time step: Euler if restart from rest
83         zraur  = 1. / rauw
84
85         ! Vertically integrated quantities
86         ! --------------------------------
87         zun(:,:) = 0.e0
88         zvn(:,:) = 0.e0
89         !
90         DO jk = 1, jpkm1             ! Vertically integrated transports (now)
91            zun(:,:) = zun(:,:) + fse3u(:,:,jk) * un(:,:,jk)
92            zvn(:,:) = zvn(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
93         END DO
94
95         ! Horizontal divergence of barotropic transports
96         !--------------------------------------------------
97         zhdiv(:,:) = 0.e0
98         DO jj = 2, jpjm1
99            DO ji = 2, jpim1   ! vector opt.
100               zhdiv(ji,jj) = (  e2u(ji  ,jj  ) * zun(ji  ,jj  )     &
101                  &            - e2u(ji-1,jj  ) * zun(ji-1,jj  )     &
102                  &            + e1v(ji  ,jj  ) * zvn(ji  ,jj  )     &
103                  &            - e1v(ji  ,jj-1) * zvn(ji  ,jj-1) )   &
104                  &           / ( e1t(ji,jj) * e2t(ji,jj) )
105            END DO
106         END DO
107
108#if defined key_obc && ( defined key_dynspg_exp || defined key_dynspg_ts )
109         ! open boundaries (div must be zero behind the open boundary)
110         !  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
111         IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1)   = 0.e0    ! east
112         IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1)   = 0.e0    ! west
113         IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0    ! north
114         IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1)   = 0.e0    ! south
115#endif
116
117#if defined key_bdy
118         jgrd=1 !: tracer grid.
119         DO jb = 1, nblenrim(jgrd)
120           ji = nbi(jb,jgrd)
121           jj = nbj(jb,jgrd)
122           zhdiv(ji,jj) = 0.e0
123         END DO
124#endif
125
126         CALL lbc_lnk( zhdiv, 'T', 1. )
127
128         ! Sea surface elevation time stepping
129         ! -----------------------------------
130         zssha(:,:) = sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) ) * tmask(:,:,1)
131
132         ! Vertical velocity computed from bottom
133         ! --------------------------------------
134         DO jk = jpkm1, 1, -1
135            wn(:,:,jk) = wn(:,:,jk+1) - fse3t(:,:,jk) * hdivn(:,:,jk) &
136              &        - ( zssha(:,:) - sshb(:,:) ) * fsve3t(:,:,jk) * mut(:,:,jk) / z2dt
137         END DO
138
139      ELSE                             ! Fixed volume
140
141         ! Vertical velocity computed from bottom
142         ! --------------------------------------
143         DO jk = jpkm1, 1, -1
144            wn(:,:,jk) = wn(:,:,jk+1) - fse3t(:,:,jk) * hdivn(:,:,jk)
145         END DO
146     
147      ENDIF
148
149      IF(ln_ctl)   CALL prt_ctl(tab3d_1=wn, clinfo1=' w**2 -   : ', mask1=wn)
150
151   END SUBROUTINE wzv
152
153   !!======================================================================
154END MODULE wzvmod
Note: See TracBrowser for help on using the repository browser.