Index: NEMO/branches/2020/dev_r12472_ASINTER-05_Masson_CurrentFeedback/src/OCE/ZDF/zdfiwm.F90
===================================================================
--- NEMO/branches/2020/dev_r12472_ASINTER-05_Masson_CurrentFeedback/src/OCE/ZDF/zdfiwm.F90	(revision 12485)
+++ NEMO/branches/2020/dev_r12472_ASINTER-05_Masson_CurrentFeedback/src/OCE/ZDF/zdfiwm.F90	(revision 12495)
@@ -87,10 +87,10 @@
       !!              This is divided into three components:
       !!                 1. Bottom-intensified low-mode dissipation at critical slopes
-      !!                     zemx_iwm(z) = ( ecri_iwm / rau0 ) * EXP( -(H-z)/hcri_iwm )
+      !!                     zemx_iwm(z) = ( ecri_iwm / rho0 ) * EXP( -(H-z)/hcri_iwm )
       !!                                   / ( 1. - EXP( - H/hcri_iwm ) ) * hcri_iwm
       !!              where hcri_iwm is the characteristic length scale of the bottom 
       !!              intensification, ecri_iwm a map of available power, and H the ocean depth.
       !!                 2. Pycnocline-intensified low-mode dissipation
-      !!                     zemx_iwm(z) = ( epyc_iwm / rau0 ) * ( sqrt(rn2(z))^nn_zpyc )
+      !!                     zemx_iwm(z) = ( epyc_iwm / rho0 ) * ( sqrt(rn2(z))^nn_zpyc )
       !!                                   / SUM( sqrt(rn2(z))^nn_zpyc * e3w(z) )
       !!              where epyc_iwm is a map of available power, and nn_zpyc
@@ -98,5 +98,5 @@
       !!              energy dissipation.
       !!                 3. WKB-height dependent high mode dissipation
-      !!                     zemx_iwm(z) = ( ebot_iwm / rau0 ) * rn2(z) * EXP(-z_wkb(z)/hbot_iwm)
+      !!                     zemx_iwm(z) = ( ebot_iwm / rho0 ) * rn2(z) * EXP(-z_wkb(z)/hbot_iwm)
       !!                                   / SUM( rn2(z) * EXP(-z_wkb(z)/hbot_iwm) * e3w(z) )
       !!              where hbot_iwm is the characteristic length scale of the WKB bottom 
@@ -151,5 +151,5 @@
       DO_2D_11_11
          zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean
-         zfact(ji,jj) = rau0 * (  1._wp - EXP( -zhdep(ji,jj) / hcri_iwm(ji,jj) )  )
+         zfact(ji,jj) = rho0 * (  1._wp - EXP( -zhdep(ji,jj) / hcri_iwm(ji,jj) )  )
          IF( zfact(ji,jj) /= 0._wp )   zfact(ji,jj) = ecri_iwm(ji,jj) / zfact(ji,jj)
       END_2D
@@ -181,5 +181,5 @@
          !
          DO_2D_11_11
-            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) )
+            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rho0 * zfact(ji,jj) )
          END_2D
          !
@@ -196,5 +196,5 @@
          !
          DO_2D_11_11
-            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) )
+            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rho0 * zfact(ji,jj) )
          END_2D
          !
@@ -243,5 +243,5 @@
       !
       DO_2D_11_11
-         IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ebot_iwm(ji,jj) / ( rau0 * zfact(ji,jj) )
+         IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ebot_iwm(ji,jj) / ( rho0 * zfact(ji,jj) )
       END_2D
       !
@@ -255,5 +255,5 @@
       ! Calculate molecular kinematic viscosity
       znu_t(:,:,:) = 1.e-4_wp * (  17.91_wp - 0.53810_wp * ts(:,:,:,jp_tem,Kmm) + 0.00694_wp * ts(:,:,:,jp_tem,Kmm) * ts(:,:,:,jp_tem,Kmm)  &
-         &                                  + 0.02305_wp * ts(:,:,:,jp_sal,Kmm)  ) * tmask(:,:,:) * r1_rau0
+         &                                  + 0.02305_wp * ts(:,:,:,jp_sal,Kmm)  ) * tmask(:,:,:) * r1_rho0
       DO jk = 2, jpkm1
          znu_w(:,:,jk) = 0.5_wp * ( znu_t(:,:,jk-1) + znu_t(:,:,jk) ) * wmask(:,:,jk)
@@ -293,5 +293,5 @@
          END_3D
          CALL mpp_sum( 'zdfiwm', zztmp )
-         zztmp = rau0 * zztmp ! Global integral of rauo * Kz * N^2 = power contributing to mixing 
+         zztmp = rho0 * zztmp ! Global integral of rauo * Kz * N^2 = power contributing to mixing 
          !
          IF(lwp) THEN
@@ -337,5 +337,5 @@
                                     !* output useful diagnostics: Kz*N^2 , 
 !!gm Kz*N2 should take into account the ratio avs/avt if it is used.... (see diaar5)
-                                    !  vertical integral of rau0 * Kz * N^2 , energy density (zemx_iwm)
+                                    !  vertical integral of rho0 * Kz * N^2 , energy density (zemx_iwm)
       IF( iom_use("bflx_iwm") .OR. iom_use("pcmap_iwm") ) THEN
          ALLOCATE( z2d(jpi,jpj) , z3d(jpi,jpj,jpk) )
@@ -345,5 +345,5 @@
             z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)
          END DO
-         z2d(:,:) = rau0 * z2d(:,:)
+         z2d(:,:) = rho0 * z2d(:,:)
          CALL iom_put( "bflx_iwm", z3d )
          CALL iom_put( "pcmap_iwm", z2d )
