Index: branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90
===================================================================
--- branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90	(revision 4668)
+++ branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90	(revision 5086)
@@ -6,5 +6,6 @@
    !! History :  OPA  ! 2000-08  (G. Madec)  double diffusive mixing
    !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module
-   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
+   !!            3.3  ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
+   !!            3.6  ! 2013-04  (G. Madec, F. Roquet) zrau compute locally using interpolation of alpha & beta
    !!----------------------------------------------------------------------
 #if defined key_zdfddm   ||   defined key_esopa
@@ -18,4 +19,6 @@
    USE dom_oce         ! ocean space and time domain variables 
    USE zdf_oce         ! ocean vertical physics variables
+   USE eosbn2         ! equation of state
+   !
    USE in_out_manager  ! I/O manager
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
@@ -34,15 +37,15 @@
    LOGICAL , PUBLIC, PARAMETER ::   lk_zdfddm = .TRUE.  !: double diffusive mixing flag
 
-   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avs    !: salinity vertical diffusivity coeff. at w-point
-   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   rrau   !: heat/salt buoyancy flux ratio
-
-   !                      !!* Namelist namzdf_ddm : double diffusive mixing *
-   REAL(wp) ::   rn_avts   ! maximum value of avs for salt fingering
-   REAL(wp) ::   rn_hsbfr  ! heat/salt buoyancy flux ratio
+   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avs   !: salinity vertical diffusivity coeff. at w-point
+
+   !                       !!* Namelist namzdf_ddm : double diffusive mixing *
+   REAL(wp) ::   rn_avts    ! maximum value of avs for salt fingering
+   REAL(wp) ::   rn_hsbfr   ! heat/salt buoyancy flux ratio
 
    !! * Substitutions
+#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
+   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -54,6 +57,5 @@
       !!                ***  ROUTINE zdf_ddm_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( avs(jpi,jpj,jpk), rrau(jpi,jpj,jpk), STAT= zdf_ddm_alloc )
-      !
+      ALLOCATE( avs(jpi,jpj,jpk) , STAT= zdf_ddm_alloc )
       IF( lk_mpp             )   CALL mpp_sum ( zdf_ddm_alloc )
       IF( zdf_ddm_alloc /= 0 )   CALL ctl_warn('zdf_ddm_alloc: failed to allocate arrays')
@@ -71,16 +73,14 @@
       !!      diffusive mixing (i.e. salt fingering and diffusive layering)
       !!      following Merryfield et al. (1999). The rate of double diffusive 
-      !!      mixing depend on the buoyancy ratio: Rrau=alpha/beta dk[T]/dk[S]
-      !!      which is computed in rn2.F
+      !!      mixing depend on the buoyancy ratio (R=alpha/beta dk[T]/dk[S]):
       !!         * salt fingering (Schmitt 1981):
-      !!      for Rrau > 1 and rn2 > 0 : zavfs = rn_avts / ( 1 + (Rrau/rn_hsbfr)^6 )
-      !!      for Rrau > 1 and rn2 > 0 : zavfs = O
-      !!      otherwise                : zavft = 0.7 zavs / Rrau
+      !!      for R > 1 and rn2 > 0 : zavfs = rn_avts / ( 1 + (R/rn_hsbfr)^6 )
+      !!      for R > 1 and rn2 > 0 : zavfs = O
+      !!      otherwise                : zavft = 0.7 zavs / R
       !!         * diffusive layering (Federov 1988):
-      !!      for 0< Rrau < 1 and rn2 > 0 : zavdt = 1.3635e-6  
-      !!                                 * exp( 4.6 exp(-0.54 (1/Rrau-1) ) )
+      !!      for 0< R < 1 and N^2 > 0 : zavdt = 1.3635e-6 * exp( 4.6 exp(-0.54 (1/R-1) ) )
       !!      otherwise                   : zavdt = 0 
-      !!      for .5 < Rrau < 1 and rn2 > 0 : zavds = zavdt (1.885 Rrau -0.85)
-      !!      for  0 < Rrau <.5 and rn2 > 0 : zavds = zavdt 0.15 Rrau      
+      !!      for .5 < R < 1 and N^2 > 0 : zavds = zavdt (1.885 R -0.85)
+      !!      for  0 < R <.5 and N^2 > 0 : zavds = zavdt 0.15 R      
       !!      otherwise                     : zavds = 0 
       !!         * update the eddy diffusivity:
@@ -96,14 +96,16 @@
       !
       INTEGER  ::   ji, jj , jk     ! dummy loop indices
-      REAL(wp) ::   zinr, zrr       ! temporary scalars
-      REAL(wp) ::   zavft, zavfs    !    -         -
-      REAL(wp) ::   zavdt, zavds    !    -         -
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zmsks, zmskf, zmskd1, zmskd2, zmskd3
+      REAL(wp) ::   zaw, zbw, zrw   ! local scalars
+      REAL(wp) ::   zdt, zds
+      REAL(wp) ::   zinr, zrr       !   -      -
+      REAL(wp) ::   zavft, zavfs    !   -      -
+      REAL(wp) ::   zavdt, zavds    !   -      -
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('zdf_ddm')
       !
-      CALL wrk_alloc( jpi,jpj, zmsks, zmskf, zmskd1, zmskd2, zmskd3 )
-
+      CALL wrk_alloc( jpi,jpj, zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3 )
+      !
       !                                                ! ===============
       DO jk = 2, jpkm1                                 ! Horizontal slab
@@ -111,5 +113,20 @@
          ! Define the mask 
          ! ---------------
-         rrau(:,:,jk) = MAX( 1.e-20, rrau(:,:,jk) )         ! only retains positive value of rrau
+         DO jj = 1, jpj                                ! R=zrau = (alpha / beta) (dk[t] / dk[s])
+            DO ji = 1, jpi
+               zrw =   ( fsdepw(ji,jj,jk  ) - fsdept(ji,jj,jk) )   &
+                  &  / ( fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk) ) 
+               !
+               zaw = (  rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem) * zrw  )  &
+                   &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
+               zbw = (  rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal) * zrw  )  &
+                   &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
+               !
+               zdt = zaw * ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) )
+               zds = zbw * ( tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ) 
+               IF( ABS( zds) <= 1.e-20_wp )   zds = 1.e-20_wp
+               zrau(ji,jj) = MAX(  1.e-20, zdt / zds  )    ! only retains positive value of zrau
+            END DO
+         END DO
 
          DO jj = 1, jpj                                     ! indicators:
@@ -119,20 +136,20 @@
                ELSE                                       ;   zmsks(ji,jj) = 1._wp
                ENDIF
-               ! salt fingering indicator: msksf=1 if rrau>1; 0 elsewhere            
-               IF( rrau(ji,jj,jk) <= 1.          ) THEN   ;   zmskf(ji,jj) = 0._wp
+               ! salt fingering indicator: msksf=1 if R>1; 0 elsewhere            
+               IF( zrau(ji,jj) <= 1.             ) THEN   ;   zmskf(ji,jj) = 0._wp
                ELSE                                       ;   zmskf(ji,jj) = 1._wp
                ENDIF
                ! diffusive layering indicators: 
-               !     ! mskdl1=1 if 0<rrau<1; 0 elsewhere
-               IF( rrau(ji,jj,jk) >= 1.          ) THEN   ;   zmskd1(ji,jj) = 0._wp
+               !     ! mskdl1=1 if 0< R <1; 0 elsewhere
+               IF( zrau(ji,jj) >= 1.             ) THEN   ;   zmskd1(ji,jj) = 0._wp
                ELSE                                       ;   zmskd1(ji,jj) = 1._wp
                ENDIF
-               !     ! mskdl2=1 if 0<rrau<0.5; 0 elsewhere
-               IF( rrau(ji,jj,jk) >= 0.5         ) THEN   ;   zmskd2(ji,jj) = 0._wp
+               !     ! mskdl2=1 if 0< R <0.5; 0 elsewhere
+               IF( zrau(ji,jj) >= 0.5            ) THEN   ;   zmskd2(ji,jj) = 0._wp
                ELSE                                       ;   zmskd2(ji,jj) = 1._wp
                ENDIF
-               !   mskdl3=1 if 0.5<rrau<1; 0 elsewhere
-               IF( rrau(ji,jj,jk) <= 0.5 .OR. rrau(ji,jj,jk) >= 1. ) THEN   ;   zmskd3(ji,jj) = 0._wp
-               ELSE                                                         ;   zmskd3(ji,jj) = 1._wp
+               !   mskdl3=1 if 0.5< R <1; 0 elsewhere
+               IF( zrau(ji,jj) <= 0.5 .OR. zrau(ji,jj) >= 1. ) THEN   ;   zmskd3(ji,jj) = 0._wp
+               ELSE                                                   ;   zmskd3(ji,jj) = 1._wp
                ENDIF
             END DO
@@ -149,7 +166,7 @@
 !CDIR NOVERRCHK
             DO ji = 1, jpi
-               zinr = 1./rrau(ji,jj,jk)
+               zinr = 1._wp / zrau(ji,jj)
                ! salt fingering
-               zrr = rrau(ji,jj,jk)/rn_hsbfr
+               zrr = zrau(ji,jj) / rn_hsbfr
                zrr = zrr * zrr
                zavfs = rn_avts / ( 1 + zrr*zrr*zrr ) * zmsks(ji,jj) * zmskf(ji,jj)
@@ -157,6 +174,6 @@
                ! diffusive layering
                zavdt = 1.3635e-6 * EXP(  4.6 * EXP( -0.54*(zinr-1.) )  ) * zmsks(ji,jj) * zmskd1(ji,jj)
-               zavds = zavdt * zmsks(ji,jj) * (  (1.85 * rrau(ji,jj,jk) - 0.85 ) * zmskd3(ji,jj)   &
-                  &                            +  0.15 * rrau(ji,jj,jk)          * zmskd2(ji,jj)  )
+               zavds = zavdt * zmsks(ji,jj) * (  ( 1.85 * zrau(ji,jj) - 0.85 ) * zmskd3(ji,jj)   &
+                  &                             +  0.15 * zrau(ji,jj)          * zmskd2(ji,jj)  )
                ! add to the eddy viscosity coef. previously computed
                avs (ji,jj,jk) = avt(ji,jj,jk) + zavfs + zavds
@@ -196,5 +213,5 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi,jpj, zmsks, zmskf, zmskd1, zmskd2, zmskd3 )
+      CALL wrk_dealloc( jpi,jpj, zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3 )
       !
       IF( nn_timing == 1 )  CALL timing_stop('zdf_ddm')
@@ -212,6 +229,7 @@
       !!              called by zdf_ddm at the first timestep (nit000)
       !!----------------------------------------------------------------------
+      INTEGER ::   ios   ! local integer
+      !!
       NAMELIST/namzdf_ddm/ rn_avts, rn_hsbfr
-      INTEGER  ::   ios                 ! Local integer output status for namelist read
       !!----------------------------------------------------------------------
       !
