Changeset 12691 for NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/DYN/dynldf_lap_blp.F90

Timestamp:

2020-04-06T17:36:16+02:00 (3 years ago)

Author:

davestorkey

Message:

UKMO/NEMO_4.0.1_biharmonic_GM : stabilising correction terms

File:

• NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/DYN/dynldf_lap_blp.F90 (modified) (5 diffs)

NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/DYN/dynldf_lap_blp.F90

@@ -17 +17 @@
    USE ldfslp         ! iso-neutral slopes 
    USE zdf_oce        ! ocean vertical physics
+   USE phycst
    !
    USE in_out_manager ! I/O manager
@@ -211 +212 @@
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       !
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zulap, zvlap   ! laplacian at u- and v-point
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zulapdz, zvlapdz   ! -1*bhm * d/dz(del^2 u) at u- and v-point
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zulap, zvlap     ! laplacian at u- and v-point
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zulapdz, zvlapdz ! -1*bhm * d/dz(del^2 u) at u- and v-point
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zmu              ! = bhm / avm
       !!----------------------------------------------------------------------
       !
@@ -222 +224 @@
       !
   
+      ! Calculate ratio of bhm and avm for stabilising correction terms
+      ! and add to avm to be included in the vertical diffusion calculation later.
+      DO jk = 2, jpkm1 
+         DO jj = 2, jpjm1
+            DO ji = fs_2, jpim1   ! vector opt.
+               zmu(ji,jj,jk) = ( bhm(ji,jj,jk) / (avm(ji,jj,jk) + rsmall) ) * wmask(ji,jj,jk)
+               avm(ji,jj,jk) = avm(ji,jj,jk) + zmu(ji,jj,jk)
+            ENDDO
+         ENDDO
+      ENDDO
+      CALL lbc_lnk_multi( 'dyn_ldf_bgm', zmu, 'W', 1. , avm, 'W', 1. ) 
+      
+      ! Calculate (del2 u)
       CALL dyn_ldf_lap_no_ahm( kt, pub, pvb, zulap, zvlap, 1 )   
 
@@ -230 +245 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, jpim1   ! vector opt.
-                zulapdz(ji,jj,jk) = -0.5_wp*(bhm(ji+1,jj  ,jk)+bhm(ji,jj,jk))*(zulap(ji,jj,jk-1) - zulap(ji,jj,jk)) * wumask(ji,jj,jk) / e3uw_n(ji,jj,jk)
-                zvlapdz(ji,jj,jk) = -0.5_wp*(bhm(ji  ,jj+1,jk)+bhm(ji,jj,jk))*(zvlap(ji,jj,jk-1) - zvlap(ji,jj,jk)) * wvmask(ji,jj,jk) / e3vw_n(ji,jj,jk)
+                zulapdz(ji,jj,jk) = ( -0.5_wp*(bhm(ji+1,jj  ,jk)+bhm(ji,jj,jk))*(zulap(ji,jj,jk-1) - zulap(ji,jj,jk)) &
+     &                                -0.5_wp*(zmu(ji+1,jj  ,jk)+zmu(ji,jj,jk))*(pub  (ji,jj,jk-1) - pub  (ji,jj,jk)) ) * wumask(ji,jj,jk) / e3uw_n(ji,jj,jk) 
+                zvlapdz(ji,jj,jk) = ( -0.5_wp*(bhm(ji  ,jj+1,jk)+bhm(ji,jj,jk))*(zvlap(ji,jj,jk-1) - zvlap(ji,jj,jk)) &
+     &                                -0.5_wp*(zmu(ji  ,jj+1,jk)+zmu(ji,jj,jk))*(pvb  (ji,jj,jk-1) - pvb  (ji,jj,jk)) ) * wumask(ji,jj,jk) / e3uw_n(ji,jj,jk) 
             ENDDO
          ENDDO
       ENDDO
-
 
 !CW: set boundary conditions: d/dz(del^2 u) = 0 at top and bottom, so that eddy-induced velocity, w*=0
@@ -252 +268 @@
         ENDDO
 
-
-
-
 !!    calculate d/dz(-bhm * d/dz(del^2 u))
                                                        ! ===============