Changeset 592 for trunk/NEMO/OPA_SRC/DYN/wzvmod.F90

Timestamp:

2007-02-09T10:15:25+01:00 (17 years ago)

Author:

opalod

Message:

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

File:

• trunk/NEMO/OPA_SRC/DYN/wzvmod.F90 (modified) (7 diffs)

trunk/NEMO/OPA_SRC/DYN/wzvmod.F90

@@ -4 +4 @@
    !! Ocean diagnostic variable : vertical velocity
    !!==============================================================================
-
+   !! History :   5.0  !  90-10  (C. Levy, G. Madec)  Original code
+   !!             7.0  !  96-01  (G. Madec)  Statement function for e3
+   !!             8.5  !  02-07  (G. Madec)  Free form, F90
    !!----------------------------------------------------------------------
    !!   wzv        : Compute the vertical velocity
@@ -14 +16 @@
    USE prtctl          ! Print control
 
+   USE domvvl          ! Variable volume
+   USE phycst
+   USE ocesbc          ! ocean surface boundary condition
+   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
+
    IMPLICIT NONE
    PRIVATE
 
    !! * Routine accessibility
-   PUBLIC wzv       ! routine called by step.F90 and inidtr.F90
+   PUBLIC wzv          ! routine called by step.F90 and inidtr.F90
 
    !! * Substitutions
 #  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !!  OPA 9.0 , LOCEAN-IPSL (2005) 
+   !! $Header$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -40 +51 @@
       !!     velocity is computed by integrating the horizontal divergence 
       !!     from the bottom to the surface.
-      !!       The boundary conditions are w=0 at the bottom (no flux) and,
-      !!     in regid-lid case, w=0 at the sea surface.
+      !!     The boundary conditions are w=0 at the bottom (no flux) and,
+      !!     in rigid-lid case, w=0 at the sea surface.
       !!
       !! ** action  :    wn array : the now vertical velocity
-      !!
-      !! History :
-      !!   5.0  !  90-10  (C. Levy, G. Madec)  Original code
-      !!   7.0  !  96-01  (G. Madec)  Statement function for e3
-      !!   8.5  !  02-07  (G. Madec)  Free form, F90
       !!----------------------------------------------------------------------
       !! * Arguments
@@ -55 +61 @@
       !! * Local declarations
       INTEGER ::   jj, jk      ! dummy loop indices
-      !!----------------------------------------------------------------------
-      !!  OPA 9.0 , LOCEAN-IPSL (2005) 
-      !! $Header$ 
-      !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
       !!----------------------------------------------------------------------
 
@@ -103 +105 @@
       !!
       !! ** action  :   wn array : the now vertical velocity
-      !!
-      !! History :
-      !!   9.0  !  02-07  (G. Madec)  Vector optimization
       !!----------------------------------------------------------------------
       !! * Arguments
@@ -111 +110 @@
 
       !! * Local declarations
-      INTEGER ::   jk          ! dummy loop indices
-      !!----------------------------------------------------------------------
-      !!  OPA 8.5, LODYC-IPSL (2002)
+      INTEGER  ::           jk           ! dummy loop indices
+      !! Variable volume
+      INTEGER  ::   ji, jj               ! dummy loop indices
+      REAL(wp) ::   z2dt, zraur          ! temporary scalar
+      REAL(wp), DIMENSION (jpi,jpj) ::   zssha, zun, zvn, zhdiv
       !!----------------------------------------------------------------------
 
@@ -125 +126 @@
       ENDIF
 
-      ! Computation from the bottom
-      DO jk = jpkm1, 1, -1
-         wn(:,:,jk) = wn(:,:,jk+1) - fse3t(:,:,jk) * hdivn(:,:,jk)
-      END DO
+      IF( lk_vvl ) THEN                ! Variable volume
+         !
+         z2dt = 2. * rdt                                       ! time step: leap-frog
+         IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt       ! time step: Euler if restart from rest
+         zraur  = 1. / rauw
+
+         ! Vertically integrated quantities
+         ! --------------------------------
+         zun(:,:) = 0.e0
+         zvn(:,:) = 0.e0
+         !
+         DO jk = 1, jpkm1             ! Vertically integrated transports (now)
+            zun(:,:) = zun(:,:) + fse3u(:,:,jk) * un(:,:,jk)
+            zvn(:,:) = zvn(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
+         END DO
+
+         ! Horizontal divergence of barotropic transports
+         !--------------------------------------------------
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1   ! vector opt.
+               zhdiv(ji,jj) = (  e2u(ji  ,jj  ) * zun(ji  ,jj  )     &
+                  &            - e2u(ji-1,jj  ) * zun(ji-1,jj  )     &
+                  &            + e1v(ji  ,jj  ) * zvn(ji  ,jj  )     &
+                  &            - e1v(ji  ,jj-1) * zvn(ji  ,jj-1) )   &
+                  &           / ( e1t(ji,jj) * e2t(ji,jj) )
+            END DO
+         END DO
+
+#if defined key_obc && ( key_dynspg_exp || key_dynspg_ts )
+         ! open boundaries (div must be zero behind the open boundary)
+         !  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
+         IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1)   = 0.e0    ! east
+         IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1)   = 0.e0    ! west
+         IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0    ! north
+         IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1)   = 0.e0    ! south
+#endif
+
+         CALL lbc_lnk( zhdiv, 'T', 1. )
+
+         ! Sea surface elevation time stepping
+         ! -----------------------------------
+         zssha(:,:) = sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) ) * tmask(:,:,1)
+
+         ! Vertical velocity computed from bottom
+         ! --------------------------------------
+         DO jk = jpkm1, 1, -1
+            wn(:,:,jk) = wn(:,:,jk+1) - fse3t(:,:,jk) * hdivn(:,:,jk) &
+              &        - ( zssha(:,:) - sshb(:,:) ) * fsve3t(:,:,jk) * mut(:,:,jk) / z2dt
+         END DO
+
+      ELSE                             ! Fixed volume 
+
+         ! Vertical velocity computed from bottom
+         ! --------------------------------------
+         DO jk = jpkm1, 1, -1
+            wn(:,:,jk) = wn(:,:,jk+1) - fse3t(:,:,jk) * hdivn(:,:,jk)
+         END DO
+      
+      ENDIF 
 
       IF(ln_ctl)   CALL prt_ctl(tab3d_1=wn, clinfo1=' w**2 -   : ', mask1=wn)