Index: branches/2014/dev_CNRS_CICE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfkpp.F90
===================================================================
--- branches/2014/dev_CNRS_CICE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfkpp.F90	(revision 4926)
+++ branches/2014/dev_CNRS_CICE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfkpp.F90	(revision 4933)
@@ -26,14 +26,15 @@
    USE phycst         ! physical constants
    USE eosbn2         ! equation of state
-   USE zdfddm         ! double diffusion mixing
+   USE zdfddm         ! double diffusion mixing (avs array)
+   USE lib_mpp        ! MPP library
+   USE trd_oce        ! ocean trends definition
+   USE trdtra         ! tracers trends
+   !
    USE in_out_manager ! I/O manager
-   USE lib_mpp        ! MPP library
-   USE wrk_nemo       ! work arrays
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
    USE prtctl         ! Print control
-   USE trdmod_oce     ! ocean trends definition
-   USE trdtra         ! tracers trends
+   USE wrk_nemo       ! work arrays
    USE timing         ! Timing
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
 
    IMPLICIT NONE
@@ -246,8 +247,8 @@
       REAL(wp) ::   zdelta, zdelta2, zdzup, zdzdn, zdzh, zvath, zgat1, zdat1, zkm1m, zkm1t
 #if defined key_zdfddm
-      REAL(wp) ::   zrrau, zds, zavdds, zavddt,zinr   ! double diffusion mixing
-      REAL(wp), POINTER, DIMENSION(:,:)   ::     zdifs
+      REAL(wp) ::   zrw, zkm1s                    ! local scalars
+      REAL(wp) ::   zrrau, zdt, zds, zavdds, zavddt, zinr   ! double diffusion mixing
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   zdifs
       REAL(wp), POINTER, DIMENSION(:)     ::   za2s, za3s, zkmps
-      REAL(wp) ::                            zkm1s
       REAL(wp), POINTER, DIMENSION(:,:)   ::   zblcs
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdiffus
@@ -274,19 +275,17 @@
 #endif
 
-      zviscos(:,:,:) = 0.
-      zblcm  (:,:  ) = 0. 
-      zdiffut(:,:,:) = 0.
-      zblct  (:,:  ) = 0. 
+      zviscos(:,:,:) = 0._wp
+      zblcm  (:,:  ) = 0._wp
+      zdiffut(:,:,:) = 0._wp
+      zblct  (:,:  ) = 0._wp 
 #if defined key_zdfddm
-      zdiffus(:,:,:) = 0.
-      zblcs  (:,:  ) = 0. 
-#endif
-      ghats(:,:,:) = 0.
-     
-      zBo   (:,:) = 0.
-      zBosol(:,:) = 0.
-      zustar(:,:) = 0.
-
-
+      zdiffus(:,:,:) = 0._wp
+      zblcs  (:,:  ) = 0._wp 
+#endif
+      ghats  (:,:,:) = 0._wp
+      zBo    (:,:  ) = 0._wp
+      zBosol (:,:  ) = 0._wp
+      zustar (:,:  ) = 0._wp
+      !
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! I. Interior diffusivity and viscosity at w points ( T interfaces)
@@ -332,19 +331,28 @@
                   avt (ji,jj,jk) =  avt (ji,jj,jk) + rn_difri * zfri    
                ENDIF
+               !
 #if defined key_zdfddm 
-               avs (ji,jj,jk) =  avt (ji,jj,jk)              
+               !
                !  Double diffusion mixing ; NOT IN ROUTINE ZDFDDM.F90
-               ! ------------------------------------------------------------------
-               ! only retains positive value of rrau
-               zrrau = MAX( rrau(ji,jj,jk), epsln )
-               zds   = tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal)
-               IF( zrrau > 1. .AND. zds > 0.) THEN
-                  !
-                  ! Salt fingering case.
-                  !---------------------
-                  ! Compute interior diffusivity for double diffusive mixing of
-                  ! salinity. Upper bound "zrrau" by "Rrho0"; (Rrho0=1.9, difcoefnuf=0.001).
-                  ! After that set interior diffusivity for double diffusive mixing
-                  ! of temperature
+               ! -------------------------
+               avs (ji,jj,jk) =  avt (ji,jj,jk)   
+
+               ! R=zrau = (alpha / beta) (dk[t] / dk[s])
+               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)
+               zbw = (  rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal) * zrw  ) * tmask(ji,jj,jk)
+               !
+               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
+               zrrau = MAX(  epsln , zdt / zds  )    ! only retains positive value of zrau
+               !
+               IF( zrrau > 1. .AND. zds > 0.) THEN                        ! Salt fingering case.
+                  !                                                       !---------------------
+                  ! Compute interior diffusivity for double diffusive mixing of salinity. 
+                  ! Upper bound "zrrau" by "Rrho0"; (Rrho0=1.9, difcoefnuf=0.001).
+                  ! After that set interior diffusivity for double diffusive mixing of temperature
                   zavdds = MIN( zrrau, Rrho0 )
                   zavdds = ( zavdds - 1.0 ) / ( Rrho0 - 1.0 )
@@ -353,18 +361,14 @@
                   zavdds = difssf * zavdds 
                   zavddt = 0.7 * zavdds
-               ELSEIF( zrrau < 1. .AND. zrrau > 0. .AND. zds < 0.) THEN
                   !
-                  ! Diffusive convection case.
-                  !---------------------------
-                  ! Compute interior diffusivity for double diffusive mixing of
-                  ! temperature (Marmorino and Caldwell, 1976); 
+               ELSEIF( zrrau < 1. .AND. zrrau > 0. .AND. zds < 0.) THEN   ! Diffusive convection case.
+                  !                                                       !---------------------------
+                  ! Compute interior diffusivity for double diffusive mixing of temperature (Marmorino and Caldwell, 1976); 
                   ! Compute interior diffusivity for double diffusive mixing of salinity 
                   zinr   = 1. / zrrau
                   zavddt = 0.909 * EXP( 4.6 * EXP( -0.54* ( zinr - 1. ) ) ) 
                   zavddt = difsdc * zavddt
-                  IF( zrrau < 0.5) THEN
-                     zavdds = zavddt * 0.15 * zrrau
-                  ELSE
-                     zavdds = zavddt * (1.85 * zrrau - 0.85 ) 
+                  IF( zrrau < 0.5) THEN   ;   zavdds = zavddt * 0.15 * zrrau
+                  ELSE                    ;   zavdds = zavddt * (1.85 * zrrau - 0.85 ) 
                   ENDIF
                ELSE
@@ -385,44 +389,10 @@
       !---------------------------------------------------------------------
       DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1     
-            IF( nn_eos < 1) THEN   
-               zt     = tsn(ji,jj,1,jp_tem)
-               zs     = tsn(ji,jj,1,jp_sal) - 35.0
-               zh     = fsdept(ji,jj,1)
-               !  potential volumic mass
-               zrhos  = rhop(ji,jj,1)
-               zalbet = ( ( ( - 0.255019e-07 * zt + 0.298357e-05 ) * zt   &   ! ratio alpha/beta
-                  &                               - 0.203814e-03 ) * zt   &
-                  &                               + 0.170907e-01 ) * zt   &
-                  &   + 0.665157e-01                                      &
-                  &   +     ( - 0.678662e-05 * zs                         &
-                  &           - 0.846960e-04 * zt + 0.378110e-02 ) * zs   &
-                  &   +   ( ( - 0.302285e-13 * zh                         &
-                  &           - 0.251520e-11 * zs                         &
-                  &           + 0.512857e-12 * zt * zt           ) * zh   &
-                  &           - 0.164759e-06 * zs                         &
-                  &        +(   0.791325e-08 * zt - 0.933746e-06 ) * zt   &
-                  &                               + 0.380374e-04 ) * zh
-
-               zbeta  = ( ( -0.415613e-09 * zt + 0.555579e-07 ) * zt      &   ! beta
-                  &                            - 0.301985e-05 ) * zt      &
-                  &   + 0.785567e-03                                      &
-                  &   + (     0.515032e-08 * zs                           &
-                  &         + 0.788212e-08 * zt - 0.356603e-06 ) * zs     &
-                  &   +(  (   0.121551e-17 * zh                           &
-                  &         - 0.602281e-15 * zs                           &
-                  &         - 0.175379e-14 * zt + 0.176621e-12 ) * zh     &
-                  &                             + 0.408195e-10   * zs     &
-                  &     + ( - 0.213127e-11 * zt + 0.192867e-09 ) * zt     &
-                  &                             - 0.121555e-07 ) * zh
-
-               zthermal = zbeta * zalbet / ( rcp * zrhos + epsln )
-               zhalin   = zbeta * tsn(ji,jj,1,jp_sal) * rcs
-            ELSE
-               zrhos    = rhop(ji,jj,1) + rau0 * ( 1. - tmask(ji,jj,1) )
-               zthermal = rn_alpha / ( rcp * zrhos + epsln )
-               zhalin   = rn_beta * tsn(ji,jj,1,jp_sal) * rcs
-               zbeta    = rn_beta
-            ENDIF
+         DO ji = fs_2, fs_jpim1           
+            zrhos    = rau0 * ( 1._wp + rhd(ji,jj,1) ) * tmask(ji,jj,1)
+            zthermal = rab_n(ji,jj,1,jp_tem) / ( rcp * zrhos + epsln )
+            zbeta    = rab_n(ji,jj,1,jp_sal)
+            zhalin   = zbeta * tsn(ji,jj,1,jp_sal) * rcs
+            !
             ! Radiative surface buoyancy force
             zBosol(ji,jj) = grav * zthermal * qsr(ji,jj)
@@ -435,6 +405,6 @@
             ws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * tsn(ji,jj,1,jp_sal)                          &
                &             + sfx(ji,jj)                                     ) * rcs * tmask(ji,jj,1) 
-         ENDDO
-      ENDDO
+         END DO
+      END DO
 
       zflageos = 0.5 + SIGN( 0.5, nn_eos - 1. ) 
@@ -447,6 +417,6 @@
             ! Friction velocity (zustar), at T-point : LMD94 eq. 2
             zustar(ji,jj) = SQRT( taum(ji,jj) / ( zrhos +  epsln ) )
-         ENDDO
-      ENDDO
+         END DO
+      END DO
 
 !CDIR NOVERRCHK  
@@ -1270,6 +1240,6 @@
          ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
 !!bug gm jpttdzdf ==> jpttkpp
-         CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_zdf, ztrdt )
-         CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_zdf, ztrds )
+         CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt )
+         CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds )
          DEALLOCATE( ztrdt )      ;     DEALLOCATE( ztrds )
       ENDIF
@@ -1340,5 +1310,5 @@
          IF( l_trdtrc ) THEN         ! save the non-local tracer flux trends for diagnostic
             ztrtrd(:,:,:) = tra(:,:,:,jn) - ztrtrd(:,:,:)
-            CALL trd_tra( kt, 'TRC', jn, jptra_trd_zdf, ztrtrd(:,:,:) )
+            CALL trd_tra( kt, 'TRC', jn, jptra_zdf, ztrtrd(:,:,:) )
          ENDIF
          !
@@ -1375,4 +1345,5 @@
       !!----------------------------------------------------------------------
       INTEGER  ::   ji, jj, jk     ! dummy loop indices
+      INTEGER  ::   ios            ! local integer
 #if ! defined key_kppcustom
       INTEGER  ::   jm             ! dummy loop indices     
@@ -1382,5 +1353,4 @@
       REAL(wp) ::   zustar, zucube, zustvk, zeta, zehat   ! tempory scalars
 #endif
-      INTEGER  ::   ios            ! Local integer output status for namelist read
       REAL(wp) ::   zhbf           ! tempory scalars
       LOGICAL  ::   ll_kppcustom   ! 1st ocean level taken as surface layer
