Changeset 1708 for trunk/NEMO/OPA_SRC/DYN/dynbfr.F90

Timestamp:

2009-11-04T14:24:34+01:00 (14 years ago)

Author:

rblod

Message:

Stability for explicit bottom friction, see ticket #588

File:

• trunk/NEMO/OPA_SRC/DYN/dynbfr.F90 (modified) (3 diffs)

trunk/NEMO/OPA_SRC/DYN/dynbfr.F90

@@ -4 +4 @@
    !! Ocean dynamics :  bottom friction component of the momentum mixing trend
    !!==============================================================================
-   !! History :  9.0  !  08-11  (A. C. Coward)  Original code
+   !! History :  9.0  !  2008-11  (A. C. Coward)  Original code
    !!----------------------------------------------------------------------
 
@@ -45 +45 @@
       !! ** Action  :   (ua,va)   momentum trend increased by bottom friction trend
       !!---------------------------------------------------------------------
+      USE oce, ONLY :   ztrdu => ta   ! use ta as 3D workspace
+      USE oce, ONLY :   ztrdv => sa   ! use sa as 3D workspace
+      !!
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       !! 
-      INTEGER ::   ji, jj       ! dummy loop indexes
-      INTEGER ::   ikbu, ikbv   ! temporary integers
-      INTEGER ::   ikbum1, ikbvm1
-!! + Enforce stability
-!     REAL(wp):: frdt           ! temporary real
-!! - Enforce stability
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztrdu, ztrdv   ! 3D workspace
+      INTEGER  ::   ji, jj          ! dummy loop indexes
+      INTEGER  ::   ikbu  , ikbv    ! temporary integers
+      INTEGER  ::   ikbum1, ikbvm1  !    -          -
+      REAL(wp) ::   zinv, zbfru, zbfrv   ! temporary scalar
       !!---------------------------------------------------------------------
       !
-!! + Enforce stability
-!     frdt = 0.125/rdt
-!! - Enforce stability
-      IF( l_trddyn )   THEN                      ! temporary save of ta and sa trends
+      zinv = -1. / ( 2.*rdt )
+
+      IF( l_trddyn )   THEN                      ! temporary save of ua and va trends
          ztrdu(:,:,:) = ua(:,:,:)
          ztrdv(:,:,:) = va(:,:,:)
       ENDIF
+
 # if defined key_vectopt_loop
       DO jj = 1, 1
@@ -75 +75 @@
             ikbvm1 = MAX( ikbv-1, 1 )
             !
-!! + Enforce stability
             ! Apply stability criteria
-!            bfrua(ji,jj) = MIN(bfrua(ji,jj), fse3u(ji,jj,ikbu)*frdt)
-!            bfrva(ji,jj) = MIN(bfrva(ji,jj), fse3v(ji,jj,ikbv)*frdt)
+            zbfru = MAX( bfrua(ji,jj), fse3u(ji,jj,ikbu)*zinv )
+            zbfrv = MAX( bfrva(ji,jj), fse3v(ji,jj,ikbv)*zinv )
             !
-!! - Enforce stability
-            ua(ji,jj,ikbum1) = ua(ji,jj,ikbum1) + bfrua(ji,jj) * ub(ji,jj,ikbum1) / fse3u(ji,jj,ikbu)
-            va(ji,jj,ikbvm1) = va(ji,jj,ikbvm1) + bfrva(ji,jj) * vb(ji,jj,ikbvm1) / fse3v(ji,jj,ikbu)
+            ua(ji,jj,ikbum1) = ua(ji,jj,ikbum1) + zbfru * ub(ji,jj,ikbum1) / fse3u(ji,jj,ikbu)
+            va(ji,jj,ikbvm1) = va(ji,jj,ikbvm1) + zbfrv * vb(ji,jj,ikbvm1) / fse3v(ji,jj,ikbv)
+            !
          END DO
       END DO
+
       !
       IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics