Changeset 1708 for trunk/NEMO/OPA_SRC/DYN/dynspg_ts.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/dynspg_ts.F90 (modified) (5 diffs)

trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90

@@ -94 +94 @@
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
       INTEGER  ::   icycle           ! temporary scalar
+      INTEGER  ::   ikbu, ikbv        ! temporary scalar
 
       REAL(wp) ::   zraur, zcoef, z2dt_e, z2dt_b     ! temporary scalars
@@ -102 +103 @@
       !!
       REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv, zsshb_e
+      !!
+      REAL(wp), DIMENSION(jpi,jpj) ::   zbfru  , zbfrv     ! 2D workspace
       !!
       REAL(wp), DIMENSION(jpi,jpj) ::   zsshun_e, zsshvn_e   ! 2D workspace
@@ -119 +122 @@
          IF(lwp) WRITE(numout,*) ' Number of sub cycle in 1 time-step (2 rdt) : icycle = ',  2*nn_baro
          !
-         CALL ts_rst( nit000, 'READ' )   ! read or initialize the following fields:
-         !                               ! un_b, vn_b  
+         CALL ts_rst( nit000, 'READ' )   ! read or initialize the following fields: un_b, vn_b  
          !
          ua_e  (:,:) = un_b (:,:)
@@ -131 +133 @@
             ftne(1,:) = 0.e0   ;   ftnw(1,:) = 0.e0   ;   ftse(1,:) = 0.e0   ;   ftsw(1,:) = 0.e0
             DO jj = 2, jpj
-               DO ji = 2, jpi   ! NO vector opt.
+               DO ji = fs_2, jpi   ! vector opt.
                   ftne(ji,jj) = ( ff(ji-1,jj  ) + ff(ji  ,jj  ) + ff(ji  ,jj-1) ) / 3.
                   ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj  ) + ff(ji  ,jj  ) ) / 3.
@@ -252 +254 @@
       END DO
 
-      !                                            ! Remove barotropic contribution of bottom friction 
-                                                   ! from the barotropic transport trend
+                    
+      !                                             ! Remove barotropic contribution of bottom friction 
+      !                                             ! from the barotropic transport trend
+      zcoef = -1. / z2dt_b
+# if defined key_vectopt_loop
+      DO jj = 1, 1
+         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
+# else
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1
+# endif
+            ikbu = MIN( mbathy(ji+1,jj), mbathy(ji,jj) )
+            ikbv = MIN( mbathy(ji,jj+1), mbathy(ji,jj) )
+            !
+            ! Apply stability criteria for bottom friction
+            !RBbug for vvl and external mode we may need to
+            ! use varying fse3
+            zbfru  (ji,jj) = MAX( bfrua(ji,jj), fse3u(ji,jj,ikbu)*zcoef )
+            zbfrv  (ji,jj) = MAX( bfrva(ji,jj), fse3v(ji,jj,ikbv)*zcoef )
+         END DO
+      END DO
+
       IF( lk_vvl ) THEN
-       DO jj = 2, jpjm1
-          DO ji = fs_2, fs_jpim1   ! vector opt.
-             zua(ji,jj) = zua(ji,jj) - bfrua(ji,jj) * zub(ji,jj) / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1.e0 - umask(ji,jj,1) )
-             zva(ji,jj) = zva(ji,jj) - bfrva(ji,jj) * zvb(ji,jj) / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1.e0 - vmask(ji,jj,1) )
-          END DO
-       END DO
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * zub(ji,jj)   &
+                  &       / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1.e0 - umask(ji,jj,1) )
+               zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * zvb(ji,jj)   &
+                  &       / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1.e0 - vmask(ji,jj,1) )
+            END DO
+         END DO
       ELSE
-       DO jj = 2, jpjm1
-          DO ji = fs_2, fs_jpim1   ! vector opt.
-             zua(ji,jj) = zua(ji,jj) - bfrua(ji,jj) * zub(ji,jj) * hur(ji,jj)
-             zva(ji,jj) = zva(ji,jj) - bfrva(ji,jj) * zvb(ji,jj) * hvr(ji,jj)
-          END DO
-       END DO
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * zub(ji,jj) * hur(ji,jj)
+               zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * zvb(ji,jj) * hvr(ji,jj)
+            END DO
+         END DO
       ENDIF