Changeset 3294 for trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90

Timestamp:

2012-01-28T17:44:18+01:00 (12 years ago)

Author:

rblod

Message:

Merge of 3.4beta into the trunk

File:

• trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90 (modified) (6 diffs)

trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90

@@ -2 +2 @@
    !!======================================================================
    !!                       ***  MODULE  bdydyn  ***
-   !! Unstructured Open Boundary Cond. :   Flow relaxation scheme on velocities
+   !! Unstructured Open Boundary Cond. :   Apply boundary conditions to velocities
    !!======================================================================
    !! History :  1.0  !  2005-02  (J. Chanut, A. Sellar)  Original code
@@ -10 +10 @@
    !!            3.3  !  2010-09  (E.O'Dea) modifications for Shelf configurations 
    !!            3.3  !  2010-09  (D.Storkey) add ice boundary conditions
+   !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge
    !!----------------------------------------------------------------------
 #if defined key_bdy 
@@ -15 +16 @@
    !!   'key_bdy' :                    Unstructured Open Boundary Condition
    !!----------------------------------------------------------------------
-   !!   bdy_dyn_frs    : relaxation of velocities on unstructured open boundary
-   !!   bdy_dyn_fla    : Flather condition for barotropic solution
+   !!   bdy_dyn        : split velocities into barotropic and baroclinic parts
+   !!                    and call bdy_dyn2d and bdy_dyn3d to apply boundary
+   !!                    conditions
    !!----------------------------------------------------------------------
+   USE wrk_nemo        ! Memory Allocation
+   USE timing          ! Timing
    USE oce             ! ocean dynamics and tracers 
    USE dom_oce         ! ocean space and time domain
+   USE dynspg_oce      
    USE bdy_oce         ! ocean open boundary conditions
-   USE dynspg_oce      ! for barotropic variables
-   USE phycst          ! physical constants
+   USE bdydyn2d        ! open boundary conditions for barotropic solution
+   USE bdydyn3d        ! open boundary conditions for baroclinic velocities
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE bdytides        ! for tidal harmonic forcing at boundary
    USE in_out_manager  !
 
@@ -30 +34 @@
    PRIVATE
 
-   PUBLIC   bdy_dyn_frs   ! routine called in dynspg_flt (free surface case ONLY)
-# if defined key_dynspg_exp || defined key_dynspg_ts
-   PUBLIC   bdy_dyn_fla   ! routine called in dynspg_flt (free surface case ONLY)
-# endif
+   PUBLIC   bdy_dyn     ! routine called in dynspg_flt (if lk_dynspg_flt) or 
+                        ! dyn_nxt (if lk_dynspg_ts or lk_dynspg_exp)
 
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
@@ -42 +45 @@
 CONTAINS
 
-   SUBROUTINE bdy_dyn_frs( kt )
+   SUBROUTINE bdy_dyn( kt, dyn3d_only )
       !!----------------------------------------------------------------------
-      !!                  ***  SUBROUTINE bdy_dyn_frs  ***
+      !!                  ***  SUBROUTINE bdy_dyn  ***
       !!
-      !! ** Purpose : - Apply the Flow Relaxation Scheme for dynamic in the  
-      !!                case of unstructured open boundaries.
+      !! ** Purpose : - Wrapper routine for bdy_dyn2d and bdy_dyn3d.
       !!
-      !! References :- Engedahl H., 1995: Use of the flow relaxation scheme in 
-      !!               a three-dimensional baroclinic ocean model with realistic
-      !!               topography. Tellus, 365-382.
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt   ! Main time step counter
       !!
-      INTEGER  ::   jb, jk         ! dummy loop indices
-      INTEGER  ::   ii, ij, igrd   ! local integers
-      REAL(wp) ::   zwgt           ! boundary weight
-      !!----------------------------------------------------------------------
-      !
-      IF(ln_dyn_frs) THEN       ! If this is false, then this routine does nothing. 
-         !
-         IF( kt == nit000 ) THEN
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) 'bdy_dyn_frs : Flow Relaxation Scheme on momentum'
-            IF(lwp) WRITE(numout,*) '~~~~~~~'
-         ENDIF
-         !
-         igrd = 2                      ! Relaxation of zonal velocity
-         DO jb = 1, nblen(igrd)
-            DO jk = 1, jpkm1
-               ii   = nbi(jb,igrd)
-               ij   = nbj(jb,igrd)
-               zwgt = nbw(jb,igrd)
-               ua(ii,ij,jk) = ( ua(ii,ij,jk) * ( 1.- zwgt ) + ubdy(jb,jk) * zwgt ) * umask(ii,ij,jk)
-            END DO
-         END DO
-         !
-         igrd = 3                      ! Relaxation of meridional velocity
-         DO jb = 1, nblen(igrd)
-            DO jk = 1, jpkm1
-               ii   = nbi(jb,igrd)
-               ij   = nbj(jb,igrd)
-               zwgt = nbw(jb,igrd)
-               va(ii,ij,jk) = ( va(ii,ij,jk) * ( 1.- zwgt ) + vbdy(jb,jk) * zwgt ) * vmask(ii,ij,jk)
-            END DO
-         END DO 
-         CALL lbc_lnk( ua, 'U', -1. )   ;   CALL lbc_lnk( va, 'V', -1. )   ! Boundary points should be updated
-         !
-      ENDIF ! ln_dyn_frs
-      !
-   END SUBROUTINE bdy_dyn_frs
+      INTEGER, INTENT( in )           :: kt               ! Main time step counter
+      LOGICAL, INTENT( in ), OPTIONAL :: dyn3d_only       ! T => only update baroclinic velocities
+      !!
+      INTEGER               :: jk,ii,ij,ib,igrd     ! Loop counter
+      LOGICAL               :: ll_dyn2d, ll_dyn3d  
+      !!
 
+      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn')
 
-# if defined   key_dynspg_exp   ||   defined key_dynspg_ts
-   !!----------------------------------------------------------------------
-   !!   'key_dynspg_exp'        OR              explicit sea surface height
-   !!   'key_dynspg_ts '                  split-explicit sea surface height
-   !!----------------------------------------------------------------------
-   
-!! Option to use Flather with dynspg_flt not coded yet...
+      ll_dyn2d = .true.
+      ll_dyn3d = .true.
 
-   SUBROUTINE bdy_dyn_fla( pssh )
-      !!----------------------------------------------------------------------
-      !!                 ***  SUBROUTINE bdy_dyn_fla  ***
-      !!             
-      !!              - Apply Flather boundary conditions on normal barotropic velocities 
-      !!                (ln_dyn_fla=.true. or ln_tides=.true.)
-      !!
-      !! ** WARNINGS about FLATHER implementation:
-      !!1. According to Palma and Matano, 1998 "after ssh" is used. 
-      !!   In ROMS and POM implementations, it is "now ssh". In the current 
-      !!   implementation (tested only in the EEL-R5 conf.), both cases were unstable. 
-      !!   So I use "before ssh" in the following.
-      !!
-      !!2. We assume that the normal ssh gradient at the bdy is zero. As a matter of 
-      !!   fact, the model ssh just inside the dynamical boundary is used (the outside  
-      !!   ssh in the code is not updated).
-      !!
-      !! References:  Flather, R. A., 1976: A tidal model of the northwest European
-      !!              continental shelf. Mem. Soc. R. Sci. Liege, Ser. 6,10, 141-164.     
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pssh
+      IF( PRESENT(dyn3d_only) ) THEN
+         IF( dyn3d_only ) ll_dyn2d = .false.
+      ENDIF
 
-      INTEGER  ::   jb, igrd                         ! dummy loop indices
-      INTEGER  ::   ii, ij, iim1, iip1, ijm1, ijp1   ! 2D addresses
-      REAL(wp) ::   zcorr                            ! Flather correction
-      REAL(wp) ::   zforc                            ! temporary scalar
-      !!----------------------------------------------------------------------
+      !-------------------------------------------------------
+      ! Set pointers
+      !-------------------------------------------------------
 
-      ! ---------------------------------!
-      ! Flather boundary conditions     :!
-      ! ---------------------------------! 
-     
-      IF(ln_dyn_fla .OR. ln_tides) THEN ! If these are both false, then this routine does nothing. 
+      pssh => sshn
+      phur => hur
+      phvr => hvr
+      CALL wrk_alloc(jpi,jpj,pu2d,pv2d) 
 
-         ! Fill temporary array with ssh data (here spgu):
-         igrd = 4
-         spgu(:,:) = 0.0
-         DO jb = 1, nblenrim(igrd)
-            ii = nbi(jb,igrd)
-            ij = nbj(jb,igrd)
-            IF( ln_dyn_fla ) spgu(ii, ij) = sshbdy(jb)
-            IF( ln_tides )   spgu(ii, ij) = spgu(ii, ij) + sshtide(jb)
-         END DO
-         !
-         igrd = 5      ! Flather bc on u-velocity; 
-         !             ! remember that flagu=-1 if normal velocity direction is outward
-         !             ! I think we should rather use after ssh ?
-         DO jb = 1, nblenrim(igrd)
-            ii  = nbi(jb,igrd)
-            ij  = nbj(jb,igrd) 
-            iim1 = ii + MAX( 0, INT( flagu(jb) ) )   ! T pts i-indice inside the boundary
-            iip1 = ii - MIN( 0, INT( flagu(jb) ) )   ! T pts i-indice outside the boundary 
-            !
-            zcorr = - flagu(jb) * SQRT( grav * hur_e(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) )
-            zforc = ubtbdy(jb) + utide(jb)
-            ua_e(ii,ij) = zforc + zcorr * umask(ii,ij,1) 
-         END DO
-         !
-         igrd = 6      ! Flather bc on v-velocity
-         !             ! remember that flagv=-1 if normal velocity direction is outward
-         DO jb = 1, nblenrim(igrd)
-            ii  = nbi(jb,igrd)
-            ij  = nbj(jb,igrd) 
-            ijm1 = ij + MAX( 0, INT( flagv(jb) ) )   ! T pts j-indice inside the boundary
-            ijp1 = ij - MIN( 0, INT( flagv(jb) ) )   ! T pts j-indice outside the boundary 
-            !
-            zcorr = - flagv(jb) * SQRT( grav * hvr_e(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) )
-            zforc = vbtbdy(jb) + vtide(jb)
-            va_e(ii,ij) = zforc + zcorr * vmask(ii,ij,1)
-         END DO
-         CALL lbc_lnk( ua_e, 'U', -1. )   ! Boundary points should be updated
-         CALL lbc_lnk( va_e, 'V', -1. )   !
-         !
-      ENDIF ! ln_dyn_fla .or. ln_tides
-      !
-   END SUBROUTINE bdy_dyn_fla
-#endif
+      !-------------------------------------------------------
+      ! Split velocities into barotropic and baroclinic parts
+      !-------------------------------------------------------
+
+      pu2d(:,:) = 0.e0
+      pv2d(:,:) = 0.e0
+      DO jk = 1, jpkm1   !! Vertically integrated momentum trends
+          pu2d(:,:) = pu2d(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
+          pv2d(:,:) = pv2d(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
+      END DO
+      pu2d(:,:) = pu2d(:,:) * phur(:,:)
+      pv2d(:,:) = pv2d(:,:) * phvr(:,:)
+      DO jk = 1 , jpkm1
+         ua(:,:,jk) = ua(:,:,jk) - pu2d(:,:)
+         va(:,:,jk) = va(:,:,jk) - pv2d(:,:)
+      END DO
+
+      !-------------------------------------------------------
+      ! Apply boundary conditions to barotropic and baroclinic
+      ! parts separately
+      !-------------------------------------------------------
+
+      IF( ll_dyn2d ) CALL bdy_dyn2d( kt )
+
+      IF( ll_dyn3d ) CALL bdy_dyn3d( kt )
+
+      !-------------------------------------------------------
+      ! Recombine velocities
+      !-------------------------------------------------------
+
+      DO jk = 1 , jpkm1
+         ua(:,:,jk) = ( ua(:,:,jk) + pu2d(:,:) ) * umask(:,:,jk)
+         va(:,:,jk) = ( va(:,:,jk) + pv2d(:,:) ) * vmask(:,:,jk)
+      END DO
+
+      CALL wrk_dealloc(jpi,jpj,pu2d,pv2d) 
+
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn')
+
+   END SUBROUTINE bdy_dyn
 
 #else
@@ -185 +124 @@
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE bdy_dyn_frs( kt )      ! Empty routine
-      WRITE(*,*) 'bdy_dyn_frs: You should not have seen this print! error?', kt
-   END SUBROUTINE bdy_dyn_frs
-   SUBROUTINE bdy_dyn_fla( pssh )    ! Empty routine
-      REAL :: pssh(:,:)
-      WRITE(*,*) 'bdy_dyn_fla: You should not have seen this print! error?', pssh(1,1)
-   END SUBROUTINE bdy_dyn_fla
+   SUBROUTINE bdy_dyn( kt )      ! Empty routine
+      WRITE(*,*) 'bdy_dyn: You should not have seen this print! error?', kt
+   END SUBROUTINE bdy_dyn
 #endif