Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/cfgs/AMM12/EXPREF/namelist_cfg
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/cfgs/AMM12/EXPREF/namelist_cfg	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/cfgs/AMM12/EXPREF/namelist_cfg	(revision 12443)
@@ -303,5 +303,5 @@
    ln_dynspg_ts = .true.   ! split-explicit free surface
    ln_bt_auto   = .false.  ! Number of sub-step defined from:
-   nn_baro      = 30       ! =F : the number of sub-step in rn_Dt seconds
+   nn_e      = 30       ! =F : the number of sub-step in rn_Dt seconds
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/cfgs/SHARED/namelist_ref
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/cfgs/SHARED/namelist_ref	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/cfgs/SHARED/namelist_ref	(revision 12443)
@@ -651,5 +651,5 @@
          rn_tide_gamma = 0.7                   ! Tidal tilt factor
          ln_scal_load  = .false.               ! Use scalar approximation for
-            rn_scal_load = 0.094               !     load potential
+            rn_load = 0.094               !     load potential
          ln_read_load  = .false.               ! Or read load potential from file
             cn_tide_load = 'tide_LOAD_grid_T.nc'  ! filename for load potential
@@ -989,9 +989,9 @@
       ln_bt_av      = .true.     ! Time filtering of barotropic variables
          nn_bt_flt     = 1          ! Time filter choice  = 0 None
-         !                          !                     = 1 Boxcar over   nn_baro sub-steps
-         !                          !                     = 2 Boxcar over 2*nn_baro  "    "
+         !                          !                     = 1 Boxcar over   nn_e sub-steps
+         !                          !                     = 2 Boxcar over 2*nn_e  "    "
       ln_bt_auto    = .true.     ! Number of sub-step defined from:
          rn_bt_cmax   =  0.8        ! =T : the Maximum Courant Number allowed
-         nn_baro      = 30          ! =F : the number of sub-step in rn_Dt seconds
+         nn_e      = 30          ! =F : the number of sub-step in rn_Dt seconds
       rn_bt_alpha   = 0.         ! Temporal diffusion parameter (if ln_bt_av=F)
 /
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/icedia.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/icedia.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/icedia.F90	(revision 12443)
@@ -109,7 +109,7 @@
       ! ---------------------------!
       ! they must be kept outside an IF(iom_use) because of the call to dia_rst below
-      z_frc_volbot = r1_rau0 * glob_sum( 'icedia', -( wfx_ice(:,:) + wfx_snw(:,:) + wfx_err_sub(:,:) ) * e1e2t(:,:) ) * 1.e-9   ! freshwater flux ice/snow-ocean 
-      z_frc_voltop = r1_rau0 * glob_sum( 'icedia', -( wfx_sub(:,:) + wfx_spr(:,:) )                    * e1e2t(:,:) ) * 1.e-9   ! freshwater flux ice/snow-atm
-      z_frc_sal    = r1_rau0 * glob_sum( 'icedia', -      sfx(:,:)                                     * e1e2t(:,:) ) * 1.e-9   ! salt fluxes ice/snow-ocean
+      z_frc_volbot = r1_rho0 * glob_sum( 'icedia', -( wfx_ice(:,:) + wfx_snw(:,:) + wfx_err_sub(:,:) ) * e1e2t(:,:) ) * 1.e-9   ! freshwater flux ice/snow-ocean 
+      z_frc_voltop = r1_rho0 * glob_sum( 'icedia', -( wfx_sub(:,:) + wfx_spr(:,:) )                    * e1e2t(:,:) ) * 1.e-9   ! freshwater flux ice/snow-atm
+      z_frc_sal    = r1_rho0 * glob_sum( 'icedia', -      sfx(:,:)                                     * e1e2t(:,:) ) * 1.e-9   ! salt fluxes ice/snow-ocean
       z_frc_tembot =           glob_sum( 'icedia',  qt_oce_ai(:,:)                                     * e1e2t(:,:) ) * 1.e-20  ! heat on top of ocean (and below ice)
       z_frc_temtop =           glob_sum( 'icedia',  qt_atm_oi(:,:)                                     * e1e2t(:,:) ) * 1.e-20  ! heat on top of ice-coean
@@ -137,6 +137,6 @@
       IF(  iom_use('ibgvolume') .OR. iom_use('ibgsaltco') .OR. iom_use('ibgheatco') .OR. iom_use('ibgheatfx') ) THEN
             
-         zdiff_vol = r1_rau0 * glob_sum( 'icedia', ( rhoi*vt_i(:,:) + rhos*vt_s(:,:) - vol_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9   ! freshwater trend (km3) 
-         zdiff_sal = r1_rau0 * glob_sum( 'icedia', ( rhoi*st_i(:,:)                  - sal_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9   ! salt content trend (km3*pss)
+         zdiff_vol = r1_rho0 * glob_sum( 'icedia', ( rhoi*vt_i(:,:) + rhos*vt_s(:,:) - vol_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9   ! freshwater trend (km3) 
+         zdiff_sal = r1_rho0 * glob_sum( 'icedia', ( rhoi*st_i(:,:)                  - sal_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9   ! salt content trend (km3*pss)
          zdiff_tem =           glob_sum( 'icedia', ( et_i(:,:) + et_s(:,:)           - tem_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-20  ! heat content trend (1.e20 J)
          !                               + SUM( qevap_ice * a_i_b, dim=3 )       !! clem: I think this term should not be there (but needs a check)
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/icedyn_adv_pra.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/icedyn_adv_pra.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/icedyn_adv_pra.F90	(revision 12443)
@@ -687,5 +687,5 @@
             IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
                !
-               zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos )
+               zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rho0-rhoi) * r1_rhos )
                !
                IF( zvs_excess > 0._wp ) THEN   ! snow-ice interface deplets below the ocean surface
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/icedyn_adv_umx.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/icedyn_adv_umx.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/icedyn_adv_umx.F90	(revision 12443)
@@ -1505,5 +1505,5 @@
             IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
                !
-               zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos )
+               zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rho0-rhoi) * r1_rhos )
                !
                IF( zvs_excess > 0._wp ) THEN   ! snow-ice interface deplets below the ocean surface
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/icedyn_rhg_evp.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/icedyn_rhg_evp.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/icedyn_rhg_evp.F90	(revision 12443)
@@ -116,5 +116,5 @@
       INTEGER ::   jter         ! local integers
       !
-      REAL(wp) ::   zrhoco                                              ! rau0 * rn_cio
+      REAL(wp) ::   zrhoco                                              ! rho0 * rn_cio
       REAL(wp) ::   zdtevp, z1_dtevp                                    ! time step for subcycling
       REAL(wp) ::   ecc2, z1_ecc2                                       ! square of yield ellipse eccenticity
@@ -213,5 +213,5 @@
       ! 1) define some variables and initialize arrays
       !------------------------------------------------------------------------------!
-      zrhoco = rau0 * rn_cio 
+      zrhoco = rho0 * rn_cio 
 
       ! ecc2: square of yield ellipse eccenticrity
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/iceistate.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/iceistate.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/iceistate.F90	(revision 12443)
@@ -371,6 +371,6 @@
       IF( ln_ice_embd ) THEN            ! embedded sea-ice: deplete the initial ssh below sea-ice area
          !
-         ssh(:,:,Kmm) = ssh(:,:,Kmm) - snwice_mass(:,:) * r1_rau0
-         ssh(:,:,Kbb) = ssh(:,:,Kbb) - snwice_mass(:,:) * r1_rau0
+         ssh(:,:,Kmm) = ssh(:,:,Kmm) - snwice_mass(:,:) * r1_rho0
+         ssh(:,:,Kbb) = ssh(:,:,Kbb) - snwice_mass(:,:) * r1_rho0
          !
          IF( .NOT.ln_linssh ) THEN
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/icethd.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/icethd.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/icethd.F90	(revision 12443)
@@ -116,5 +116,5 @@
       ELSE      !  if no ice dynamics => transmit directly the atmospheric stress to the ocean
          DO_2D_00_00
-            zfric(ji,jj) = r1_rau0 * SQRT( 0.5_wp *  &
+            zfric(ji,jj) = r1_rho0 * SQRT( 0.5_wp *  &
                &                         (  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   &
                &                          + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) ) * tmask(ji,jj,1)
@@ -141,10 +141,10 @@
 
          ! --- Energy needed to bring ocean surface layer until its freezing (mostly<0 but >0 if supercooling, J.m-2) --- !
-         zqfr     = rau0 * rcp * e3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) * tmask(ji,jj,1)  ! both < 0 (t_bo < sst) and > 0 (t_bo > sst)
+         zqfr     = rho0 * rcp * e3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) * tmask(ji,jj,1)  ! both < 0 (t_bo < sst) and > 0 (t_bo > sst)
          zqfr_neg = MIN( zqfr , 0._wp )                                                                    ! only < 0
 
          ! --- Sensible ocean-to-ice heat flux (mostly>0 but <0 if supercooling, W/m2)
          zfric_u            = MAX( SQRT( zfric(ji,jj) ), zfric_umin ) 
-         qsb_ice_bot(ji,jj) = rswitch * rau0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ! W.m-2
+         qsb_ice_bot(ji,jj) = rswitch * rho0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ! W.m-2
 
          qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) )
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/icethd_dh.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/icethd_dh.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/icethd_dh.F90	(revision 12443)
@@ -76,5 +76,5 @@
       REAL(wp) ::   zgrr         ! bottom growth rate
       REAL(wp) ::   zt_i_new     ! bottom formation temperature
-      REAL(wp) ::   z1_rho       ! 1/(rhos+rau0-rhoi)
+      REAL(wp) ::   z1_rho       ! 1/(rhos+rho0-rhoi)
 
       REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean
@@ -571,8 +571,8 @@
       ! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level, 
       ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice
-      z1_rho = 1._wp / ( rhos+rau0-rhoi )
+      z1_rho = 1._wp / ( rhos+rho0-rhoi )
       DO ji = 1, npti
          !
-         dh_snowice(ji) = MAX(  0._wp , ( rhos * h_s_1d(ji) + (rhoi-rau0) * h_i_1d(ji) ) * z1_rho )
+         dh_snowice(ji) = MAX(  0._wp , ( rhos * h_s_1d(ji) + (rhoi-rho0) * h_i_1d(ji) ) * z1_rho )
 
          h_i_1d(ji)    = h_i_1d(ji) + dh_snowice(ji)
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/icethd_do.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/icethd_do.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/icethd_do.F90	(revision 12443)
@@ -141,5 +141,5 @@
          ! Physical constants
          zhicrit = 0.04                                          ! frazil ice thickness
-         ztwogp  = 2. * rau0 / ( grav * 0.3 * ( rau0 - rhoi ) )  ! reduced grav
+         ztwogp  = 2. * rho0 / ( grav * 0.3 * ( rho0 - rhoi ) )  ! reduced grav
          zsqcd   = 1.0 / SQRT( 1.3 * zcai )                      ! 1/SQRT(airdensity*drag)
          zgamafr = 0.03
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/iceupdate.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/iceupdate.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/iceupdate.F90	(revision 12443)
@@ -329,5 +329,5 @@
       ENDIF
 
-      zrhoco = rau0 * rn_cio
+      zrhoco = rho0 * rn_cio
       !
       IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !==  Ice time-step only  ==!   (i.e. surface module time-step)
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/icevar.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/icevar.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/icevar.F90	(revision 12443)
@@ -717,5 +717,5 @@
       !! ** Purpose :  compute the equivalent ssh in lead when sea ice is embedded
       !!
-      !! ** Method  :  ssh_lead = ssh + (Mice + Msnow) / rau0
+      !! ** Method  :  ssh_lead = ssh + (Mice + Msnow) / rho0
       !!
       !! ** Reference : Jean-Michel Campin, John Marshall, David Ferreira,
@@ -747,5 +747,5 @@
          zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
          !
-         zsnwiceload(:,:) = ( zintn * psnwice_mass(:,:) + zintb * psnwice_mass_b(:,:) ) * r1_rau0
+         zsnwiceload(:,:) = ( zintn * psnwice_mass(:,:) + zintb * psnwice_mass_b(:,:) ) * r1_rho0
          !
       ELSE
@@ -937,5 +937,5 @@
                ! In case snow load is in excess that would lead to transformation from snow to ice
                ! Then, transfer the snow excess into the ice (different from icethd_dh)
-               zdh = MAX( 0._wp, ( rhos * ph_s(ji,jl) + ( rhoi - rau0 ) * ph_i(ji,jl) ) * r1_rau0 ) 
+               zdh = MAX( 0._wp, ( rhos * ph_s(ji,jl) + ( rhoi - rho0 ) * ph_i(ji,jl) ) * r1_rho0 ) 
                ! recompute h_i, h_s avoiding out of bounds values
                ph_i(ji,jl) = MIN( hi_max(jl), ph_i(ji,jl) + zdh )
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/icewri.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/icewri.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/ICE/icewri.F90	(revision 12443)
@@ -87,5 +87,5 @@
       ! Standard outputs
       !-----------------
-      zrho1 = ( rau0 - rhoi ) * r1_rau0 ; zrho2 = rhos * r1_rau0
+      zrho1 = ( rho0 - rhoi ) * r1_rho0 ; zrho2 = rhos * r1_rho0
       ! masks
       CALL iom_put( 'icemask'  , zmsk00 )   ! ice mask 0%
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/NST/agrif_oce_update.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/NST/agrif_oce_update.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/NST/agrif_oce_update.F90	(revision 12443)
@@ -361,5 +361,5 @@
                            ztnu = tabres_child(ji,jj,jk,jn) * e3t(ji,jj,jk,Kmm_a)
                            ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a)
-                           ts(ji,jj,jk,jn,Kbb_a) = ( ztb + atfp * ( ztnu - ztno) )  & 
+                           ts(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) )  & 
                                      &        * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a)
                         ENDIF
@@ -432,5 +432,5 @@
                            ztnu = tabres(ji,jj,jk,jn)
                            ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a)
-                           ts(ji,jj,jk,jn,Kbb_a) = ( ztb + atfp * ( ztnu - ztno) )  & 
+                           ts(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) )  & 
                                      &        * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a)
                         ENDIF
@@ -555,5 +555,5 @@
                      zuno = uu(ji,jj,jk,Kmm_a) * e3u(ji,jj,jk,Krhs_a)
                      zunu = tabres_child(ji,jj,jk) * e3u(ji,jj,jk,Kmm_a)
-                     uu(ji,jj,jk,Kbb_a) = ( zub + atfp * ( zunu - zuno) ) &      
+                     uu(ji,jj,jk,Kbb_a) = ( zub + rn_atfp * ( zunu - zuno) ) &      
                                     & * umask(ji,jj,jk) / e3u(ji,jj,jk,Kbb_a)
                   ENDIF
@@ -601,5 +601,5 @@
                      zuno = uu(ji,jj,jk,Kmm_a) * e3u(ji,jj,jk,Krhs_a)
                      zunu = tabres(ji,jj,jk,1)
-                     uu(ji,jj,jk,Kbb_a) = ( zub + atfp * ( zunu - zuno) ) &      
+                     uu(ji,jj,jk,Kbb_a) = ( zub + rn_atfp * ( zunu - zuno) ) &      
                                     & * umask(ji,jj,jk) / e3u(ji,jj,jk,Kbb_a)
                   ENDIF
@@ -755,5 +755,5 @@
                      zvno = vv(ji,jj,jk,Kmm_a) * e3v(ji,jj,jk,Krhs_a)
                      zvnu = tabres_child(ji,jj,jk) * e3v(ji,jj,jk,Kmm_a)
-                     vv(ji,jj,jk,Kbb_a) = ( zvb + atfp * ( zvnu - zvno) ) &      
+                     vv(ji,jj,jk,Kbb_a) = ( zvb + rn_atfp * ( zvnu - zvno) ) &      
                                     & * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kbb_a)
                   ENDIF
@@ -805,5 +805,5 @@
                      zvno = vv(ji,jj,jk,Kmm_a) * e3v(ji,jj,jk,Krhs_a)
                      zvnu = tabres(ji,jj,jk,1)
-                     vv(ji,jj,jk,Kbb_a) = ( zvb + atfp * ( zvnu - zvno) ) &      
+                     vv(ji,jj,jk,Kbb_a) = ( zvb + rn_atfp * ( zvnu - zvno) ) &      
                                     & * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kbb_a)
                   ENDIF
@@ -909,5 +909,5 @@
                   IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part
                      zcorr = (tabres(ji,jj) - uu_b(ji,jj,Kmm_a) * hu(ji,jj,Krhs_a)) * r1_hu(ji,jj,Kbb_a)
-                     uu_b(ji,jj,Kbb_a) = uu_b(ji,jj,Kbb_a) + atfp * zcorr * umask(ji,jj,1)
+                     uu_b(ji,jj,Kbb_a) = uu_b(ji,jj,Kbb_a) + rn_atfp * zcorr * umask(ji,jj,1)
                   END IF
                ENDIF    
@@ -975,5 +975,5 @@
                   IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part
                      zcorr = (tabres(ji,jj) - vv_b(ji,jj,Kmm_a) * hv(ji,jj,Krhs_a)) * r1_hv(ji,jj,Kbb_a)
-                     vv_b(ji,jj,Kbb_a) = vv_b(ji,jj,Kbb_a) + atfp * zcorr * vmask(ji,jj,1)
+                     vv_b(ji,jj,Kbb_a) = vv_b(ji,jj,Kbb_a) + rn_atfp * zcorr * vmask(ji,jj,1)
                   END IF
                ENDIF              
@@ -1025,5 +1025,5 @@
                DO ji=i1,i2
                   ssh(ji,jj,Kbb_a) =   ssh(ji,jj,Kbb_a) &
-                        & + atfp * ( tabres(ji,jj) - ssh(ji,jj,Kmm_a) ) * tmask(ji,jj,1)
+                        & + rn_atfp * ( tabres(ji,jj) - ssh(ji,jj,Kmm_a) ) * tmask(ji,jj,1)
                END DO
             END DO
@@ -1119,5 +1119,5 @@
                zcor = rn_Dt * r1_e1e2t(i1  ,jj) * e2u(i1,jj) * (ub2_b(i1,jj)-tabres(i1,jj)) 
                ssh(i1  ,jj,Kmm_a) = ssh(i1  ,jj,Kmm_a) + zcor
-               IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) ssh(i1  ,jj,Kbb_a) = ssh(i1  ,jj,Kbb_a) + atfp * zcor
+               IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) ssh(i1  ,jj,Kbb_a) = ssh(i1  ,jj,Kbb_a) + rn_atfp * zcor
             END DO
          ENDIF
@@ -1126,5 +1126,5 @@
                zcor = - rn_Dt * r1_e1e2t(i2+1,jj) * e2u(i2,jj) * (ub2_b(i2,jj)-tabres(i2,jj))
                ssh(i2+1,jj,Kmm_a) = ssh(i2+1,jj,Kmm_a) + zcor
-               IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) ssh(i2+1,jj,Kbb_a) = ssh(i2+1,jj,Kbb_a) + atfp * zcor
+               IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) ssh(i2+1,jj,Kbb_a) = ssh(i2+1,jj,Kbb_a) + rn_atfp * zcor
             END DO
          ENDIF
@@ -1207,5 +1207,5 @@
                zcor = rn_Dt * r1_e1e2t(ji,j1  ) * e1v(ji,j1  ) * (vb2_b(ji,j1)-tabres(ji,j1))
                ssh(ji,j1  ,Kmm_a) = ssh(ji,j1  ,Kmm_a) + zcor
-               IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) ssh(ji,j1  ,Kbb_a) = ssh(ji,j1,Kbb_a) + atfp * zcor
+               IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) ssh(ji,j1  ,Kbb_a) = ssh(ji,j1,Kbb_a) + rn_atfp * zcor
             END DO
          ENDIF
@@ -1214,5 +1214,5 @@
                zcor = - rn_Dt * r1_e1e2t(ji,j2+1) * e1v(ji,j2  ) * (vb2_b(ji,j2)-tabres(ji,j2))
                ssh(ji,j2+1,Kmm_a) = ssh(ji,j2+1,Kmm_a) + zcor
-               IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) ssh(ji,j2+1,Kbb_a) = ssh(ji,j2+1,Kbb_a) + atfp * zcor
+               IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) ssh(ji,j2+1,Kbb_a) = ssh(ji,j2+1,Kbb_a) + rn_atfp * zcor
             END DO
          ENDIF
@@ -1364,5 +1364,5 @@
                   DO ji=i1,i2
                      e3t(ji,jj,jk,Kbb_a) =  e3t(ji,jj,jk,Kbb_a) &
-                           & + atfp * ( ptab(ji,jj,jk) - e3t(ji,jj,jk,Kmm_a) )
+                           & + rn_atfp * ( ptab(ji,jj,jk) - e3t(ji,jj,jk,Kmm_a) )
                   END DO
                END DO
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/NST/agrif_top_update.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/NST/agrif_top_update.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/NST/agrif_top_update.F90	(revision 12443)
@@ -135,5 +135,5 @@
                            ztnu = tabres_child(ji,jj,jk,jn) * e3t(ji,jj,jk,Kmm_a)
                            ztno = tr(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a)
-                           tr(ji,jj,jk,jn,Kbb_a) = ( ztb + atfp * ( ztnu - ztno) )  & 
+                           tr(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) )  & 
                                      &        * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a)
                         ENDIF
@@ -209,5 +209,5 @@
                            ztnu = tabres(ji,jj,jk,jn)
                            ztno = tr(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a)
-                           tr(ji,jj,jk,jn,Kbb_a) = ( ztb + atfp * ( ztnu - ztno) )  & 
+                           tr(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) )  & 
                                      &        * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a)
                         ENDIF
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/BDY/bdyice.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/BDY/bdyice.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/BDY/bdyice.F90	(revision 12443)
@@ -179,5 +179,5 @@
 
             ! Then, a) transfer the snow excess into the ice (different from icethd_dh)
-            zdh = MAX( 0._wp, ( rhos * h_s(ji,jj,jl) + ( rhoi - rau0 ) * h_i(ji,jj,jl) ) * r1_rau0 )
+            zdh = MAX( 0._wp, ( rhos * h_s(ji,jj,jl) + ( rhoi - rho0 ) * h_i(ji,jj,jl) ) * r1_rho0 )
             ! Or, b) transfer all the snow into ice (if incoming ice is likely to melt as it comes into a warmer environment)
             !zdh = MAX( 0._wp, h_s(ji,jj,jl) * rhos / rhoi )
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/BDY/bdytides.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/BDY/bdytides.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/BDY/bdytides.F90	(revision 12443)
@@ -297,5 +297,5 @@
       ! Absolute time from model initialization:   
       IF( PRESENT(kit) ) THEN  
-         z_arg = ( REAL(kt, wp) + ( REAL(kit, wp) + zt_offset - 1. ) / REAL(nn_baro, wp) ) * rn_Dt
+         z_arg = ( REAL(kt, wp) + ( REAL(kit, wp) + zt_offset - 1. ) / REAL(nn_e, wp) ) * rn_Dt
       ELSE                              
          z_arg = ( REAL(kt, wp) + zt_offset ) * rn_Dt
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/BDY/bdyvol.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/BDY/bdyvol.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/BDY/bdyvol.F90	(revision 12443)
@@ -77,5 +77,5 @@
       ! Calculate the cumulate surface Flux z_cflxemp (m3/s) over all the domain
       ! -----------------------------------------------------------------------
-      IF ( kc == 1 ) z_cflxemp = glob_sum( 'bdyvol', ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) ) * bdytmask(:,:) * e1e2t(:,:)  ) / rau0
+      IF ( kc == 1 ) z_cflxemp = glob_sum( 'bdyvol', ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) ) * bdytmask(:,:) * e1e2t(:,:)  ) / rho0
 
       ! Compute bdy surface each cycle if non linear free surface
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DIA/diaar5.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DIA/diaar5.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DIA/diaar5.F90	(revision 12443)
@@ -103,5 +103,5 @@
          END DO
          CALL iom_put( 'volcello'  , zrhd(:,:,:)  )  ! WARNING not consistent with CMIP DR where volcello is at ca. 2000
-         CALL iom_put( 'masscello' , rau0 * e3t(:,:,:,Kmm) * tmask(:,:,:) )  ! ocean mass
+         CALL iom_put( 'masscello' , rho0 * e3t(:,:,:,Kmm) * tmask(:,:,:) )  ! ocean mass
       ENDIF 
       !
@@ -181,5 +181,5 @@
          CALL iom_put( 'sshsteric', zssh_steric )
          !                                         ! ocean bottom pressure
-         zztmp = rau0 * grav * 1.e-4_wp               ! recover pressure from pressure anomaly and cover to dbar = 1.e4 Pa
+         zztmp = rho0 * grav * 1.e-4_wp               ! recover pressure from pressure anomaly and cover to dbar = 1.e4 Pa
          zbotpres(:,:) = zztmp * ( zbotpres(:,:) + ssh(:,:,Kmm) + thick0(:,:) )
          CALL iom_put( 'botpres', zbotpres )
@@ -213,5 +213,5 @@
          ztemp = glob_sum( 'diaar5', ztsn(:,:,1,jp_tem) )
          zsal  = glob_sum( 'diaar5', ztsn(:,:,1,jp_sal) )
-         zmass = rau0 * ( zarho + zvol )      
+         zmass = rho0 * ( zarho + zvol )      
          !
          CALL iom_put( 'masstot', zmass )
@@ -251,5 +251,5 @@
                z2d(:,:) = 0._wp
                DO_3D_11_11( 1, jpkm1 )
-                  z2d(ji,jj) = z2d(ji,jj) + rau0 * e3t(ji,jj,jk,Kmm) *  ztpot(ji,jj,jk)
+                  z2d(ji,jj) = z2d(ji,jj) + rho0 * e3t(ji,jj,jk,Kmm) *  ztpot(ji,jj,jk)
                END_3D
                CALL iom_put( 'tosmint_pot', z2d ) 
@@ -285,5 +285,5 @@
           ELSE
             DO_3D_11_11( 1, jpk )
-               zpe(ji,jj) = zpe(ji,jj) + avt(ji,jj,jk) * MIN(0._wp,rn2(ji,jj,jk)) * rau0 * e3w(ji,jj,jk,Kmm)
+               zpe(ji,jj) = zpe(ji,jj) + avt(ji,jj,jk) * MIN(0._wp,rn2(ji,jj,jk)) * rho0 * e3w(ji,jj,jk,Kmm)
             END_3D
          ENDIF
@@ -325,10 +325,10 @@
        CALL lbc_lnk( 'diaar5', z2d, 'U', -1. )
        IF( cptr == 'adv' ) THEN
-          IF( ktra == jp_tem ) CALL iom_put( 'uadv_heattr' , rau0_rcp * z2d )  ! advective heat transport in i-direction
-          IF( ktra == jp_sal ) CALL iom_put( 'uadv_salttr' , rau0     * z2d )  ! advective salt transport in i-direction
+          IF( ktra == jp_tem ) CALL iom_put( 'uadv_heattr' , rho0_rcp * z2d )  ! advective heat transport in i-direction
+          IF( ktra == jp_sal ) CALL iom_put( 'uadv_salttr' , rho0     * z2d )  ! advective salt transport in i-direction
        ENDIF
        IF( cptr == 'ldf' ) THEN
-          IF( ktra == jp_tem ) CALL iom_put( 'udiff_heattr' , rau0_rcp * z2d ) ! diffusive heat transport in i-direction
-          IF( ktra == jp_sal ) CALL iom_put( 'udiff_salttr' , rau0     * z2d ) ! diffusive salt transport in i-direction
+          IF( ktra == jp_tem ) CALL iom_put( 'udiff_heattr' , rho0_rcp * z2d ) ! diffusive heat transport in i-direction
+          IF( ktra == jp_sal ) CALL iom_put( 'udiff_salttr' , rho0     * z2d ) ! diffusive salt transport in i-direction
        ENDIF
        !
@@ -339,10 +339,10 @@
        CALL lbc_lnk( 'diaar5', z2d, 'V', -1. )
        IF( cptr == 'adv' ) THEN
-          IF( ktra == jp_tem ) CALL iom_put( 'vadv_heattr' , rau0_rcp * z2d )  ! advective heat transport in j-direction
-          IF( ktra == jp_sal ) CALL iom_put( 'vadv_salttr' , rau0     * z2d )  ! advective salt transport in j-direction
+          IF( ktra == jp_tem ) CALL iom_put( 'vadv_heattr' , rho0_rcp * z2d )  ! advective heat transport in j-direction
+          IF( ktra == jp_sal ) CALL iom_put( 'vadv_salttr' , rho0     * z2d )  ! advective salt transport in j-direction
        ENDIF
        IF( cptr == 'ldf' ) THEN
-          IF( ktra == jp_tem ) CALL iom_put( 'vdiff_heattr' , rau0_rcp * z2d ) ! diffusive heat transport in j-direction
-          IF( ktra == jp_sal ) CALL iom_put( 'vdiff_salttr' , rau0     * z2d ) ! diffusive salt transport in j-direction
+          IF( ktra == jp_tem ) CALL iom_put( 'vdiff_heattr' , rho0_rcp * z2d ) ! diffusive heat transport in j-direction
+          IF( ktra == jp_sal ) CALL iom_put( 'vdiff_salttr' , rho0     * z2d ) ! diffusive salt transport in j-direction
        ENDIF
           
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DIA/diadct.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DIA/diadct.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DIA/diadct.F90	(revision 12443)
@@ -676,5 +676,5 @@
                   zsn   = interp(Kmm,k%I,k%J,jk,'V',ts(:,:,:,jp_sal,Kmm) ) 
                   zrhop = interp(Kmm,k%I,k%J,jk,'V',rhop) 
-                  zrhoi = interp(Kmm,k%I,k%J,jk,'V',rhd*rau0+rau0) 
+                  zrhoi = interp(Kmm,k%I,k%J,jk,'V',rhd*rho0+rho0) 
                   zsshn =  0.5*( ssh(k%I,k%J,Kmm) + ssh(k%I,k%J+1,Kmm)    ) * vmask(k%I,k%J,1) 
                CASE(2,3) 
@@ -682,5 +682,5 @@
                   zsn   = interp(Kmm,k%I,k%J,jk,'U',ts(:,:,:,jp_sal,Kmm) ) 
                   zrhop = interp(Kmm,k%I,k%J,jk,'U',rhop) 
-                  zrhoi = interp(Kmm,k%I,k%J,jk,'U',rhd*rau0+rau0) 
+                  zrhoi = interp(Kmm,k%I,k%J,jk,'U',rhd*rho0+rho0) 
                   zsshn =  0.5*( ssh(k%I,k%J,Kmm) + ssh(k%I+1,k%J,Kmm)    ) * umask(k%I,k%J,1)  
                END SELECT 
@@ -849,5 +849,5 @@
                  zsn   = interp(Kmm,k%I,k%J,jk,'V',ts(:,:,:,jp_sal,Kmm) ) 
                  zrhop = interp(Kmm,k%I,k%J,jk,'V',rhop) 
-                 zrhoi = interp(Kmm,k%I,k%J,jk,'V',rhd*rau0+rau0) 
+                 zrhoi = interp(Kmm,k%I,k%J,jk,'V',rhd*rho0+rho0) 
 
               CASE(2,3) 
@@ -855,5 +855,5 @@
                  zsn   = interp(Kmm,k%I,k%J,jk,'U',ts(:,:,:,jp_sal,Kmm) ) 
                  zrhop = interp(Kmm,k%I,k%J,jk,'U',rhop) 
-                 zrhoi = interp(Kmm,k%I,k%J,jk,'U',rhd*rau0+rau0) 
+                 zrhoi = interp(Kmm,k%I,k%J,jk,'U',rhd*rho0+rho0) 
                  zsshn =  0.5*( ssh(k%I,k%J,Kmm)    + ssh(k%I+1,k%J,Kmm)    ) * umask(k%I,k%J,1)  
               END SELECT 
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DIA/diahsb.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DIA/diahsb.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DIA/diahsb.F90	(revision 12443)
@@ -91,5 +91,5 @@
       ! 1 - Trends due to forcing !
       ! ------------------------- !
-      z_frc_trd_v = r1_rau0 * glob_sum( 'diahsb', - ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) ) * surf(:,:) )   ! volume fluxes
+      z_frc_trd_v = r1_rho0 * glob_sum( 'diahsb', - ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) ) * surf(:,:) )   ! volume fluxes
       z_frc_trd_t =           glob_sum( 'diahsb', sbc_tsc(:,:,jp_tem) * surf(:,:) )                       ! heat fluxes
       z_frc_trd_s =           glob_sum( 'diahsb', sbc_tsc(:,:,jp_sal) * surf(:,:) )                       ! salt fluxes
@@ -101,5 +101,5 @@
          &                          + glob_sum( 'diahsb', ( risf_cav_tsc(:,:,jp_tem) + risf_par_tsc(:,:,jp_tem) ) * surf(:,:) )
       !                    ! Add penetrative solar radiation
-      IF( ln_traqsr )   z_frc_trd_t = z_frc_trd_t + r1_rau0_rcp * glob_sum( 'diahsb', qsr     (:,:) * surf(:,:) )
+      IF( ln_traqsr )   z_frc_trd_t = z_frc_trd_t + r1_rho0_rcp * glob_sum( 'diahsb', qsr     (:,:) * surf(:,:) )
       !                    ! Add geothermal heat flux
       IF( ln_trabbc )   z_frc_trd_t = z_frc_trd_t +               glob_sum( 'diahsb', qgh_trd0(:,:) * surf(:,:) )
@@ -197,6 +197,6 @@
 
       CALL iom_put(   'bgfrcvol' , frc_v    * 1.e-9    )              ! vol - surface forcing (km3) 
-      CALL iom_put(   'bgfrctem' , frc_t    * rau0 * rcp * 1.e-20 )   ! hc  - surface forcing (1.e20 J) 
-      CALL iom_put(   'bgfrchfx' , frc_t    * rau0 * rcp /  &         ! hc  - surface forcing (W/m2) 
+      CALL iom_put(   'bgfrctem' , frc_t    * rho0 * rcp * 1.e-20 )   ! hc  - surface forcing (1.e20 J) 
+      CALL iom_put(   'bgfrchfx' , frc_t    * rho0 * rcp /  &         ! hc  - surface forcing (W/m2) 
          &                       ( surf_tot * kt * rn_Dt )        )
       CALL iom_put(   'bgfrcsal' , frc_s    * 1.e-9    )              ! sc  - surface forcing (psu*km3) 
@@ -205,6 +205,6 @@
          CALL iom_put( 'bgtemper' , zdiff_hc / zvol_tot )              ! Temperature drift     (C) 
          CALL iom_put( 'bgsaline' , zdiff_sc / zvol_tot )              ! Salinity    drift     (PSU)
-         CALL iom_put( 'bgheatco' , zdiff_hc * 1.e-20 * rau0 * rcp )   ! Heat content drift    (1.e20 J) 
-         CALL iom_put( 'bgheatfx' , zdiff_hc * rau0 * rcp /  &         ! Heat flux drift       (W/m2) 
+         CALL iom_put( 'bgheatco' , zdiff_hc * 1.e-20 * rho0 * rcp )   ! Heat content drift    (1.e20 J) 
+         CALL iom_put( 'bgheatfx' , zdiff_hc * rho0 * rcp /  &         ! Heat flux drift       (W/m2) 
             &                       ( surf_tot * kt * rn_Dt )        )
          CALL iom_put( 'bgsaltco' , zdiff_sc * 1.e-9    )              ! Salt content drift    (psu*km3)
@@ -225,6 +225,6 @@
          CALL iom_put( 'bgtemper' , zdiff_hc1 / zvol_tot)              ! Heat content drift    (C) 
          CALL iom_put( 'bgsaline' , zdiff_sc1 / zvol_tot)              ! Salt content drift    (PSU)
-         CALL iom_put( 'bgheatco' , zdiff_hc1 * 1.e-20 * rau0 * rcp )  ! Heat content drift    (1.e20 J) 
-         CALL iom_put( 'bgheatfx' , zdiff_hc1 * rau0 * rcp /  &        ! Heat flux drift       (W/m2) 
+         CALL iom_put( 'bgheatco' , zdiff_hc1 * 1.e-20 * rho0 * rcp )  ! Heat content drift    (1.e20 J) 
+         CALL iom_put( 'bgheatfx' , zdiff_hc1 * rho0 * rcp /  &        ! Heat flux drift       (W/m2) 
             &                       ( surf_tot * kt * rn_Dt )         )
          CALL iom_put( 'bgsaltco' , zdiff_sc1 * 1.e-9    )             ! Salt content drift    (psu*km3)
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DIA/diahth.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DIA/diahth.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DIA/diahth.F90	(revision 12443)
@@ -261,5 +261,5 @@
             zzdep = 300.
             CALL  dia_hth_htc( Kmm, zzdep, ts(:,:,:,jp_tem,Kmm), htc3 )
-            CALL iom_put( 'hc300', rau0_rcp * htc3 )  ! vertically integrated heat content (J/m2)
+            CALL iom_put( 'hc300', rho0_rcp * htc3 )  ! vertically integrated heat content (J/m2)
          ENDIF
          !
@@ -270,5 +270,5 @@
             zzdep = 700.
             CALL  dia_hth_htc( Kmm, zzdep, ts(:,:,:,jp_tem,Kmm), htc7 )
-            CALL iom_put( 'hc700', rau0_rcp * htc7 )  ! vertically integrated heat content (J/m2)
+            CALL iom_put( 'hc700', rho0_rcp * htc7 )  ! vertically integrated heat content (J/m2)
   
          ENDIF
@@ -280,5 +280,5 @@
             zzdep = 2000.
             CALL  dia_hth_htc( Kmm, zzdep, ts(:,:,:,jp_tem,Kmm), htc20 )
-            CALL iom_put( 'hc2000', rau0_rcp * htc20 )  ! vertically integrated heat content (J/m2)  
+            CALL iom_put( 'hc2000', rho0_rcp * htc20 )  ! vertically integrated heat content (J/m2)  
          ENDIF
          !
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DIA/diaptr.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DIA/diaptr.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DIA/diaptr.F90	(revision 12443)
@@ -50,6 +50,6 @@
 
    REAL(wp) ::   rc_sv    = 1.e-6_wp   ! conversion from m3/s to Sverdrup
-   REAL(wp) ::   rc_pwatt = 1.e-15_wp  ! conversion from W    to PW (further x rau0 x Cp)
-   REAL(wp) ::   rc_ggram = 1.e-9_wp   ! conversion from g    to Gg  (further x rau0)
+   REAL(wp) ::   rc_pwatt = 1.e-15_wp  ! conversion from W    to PW (further x rho0 x Cp)
+   REAL(wp) ::   rc_ggram = 1.e-9_wp   ! conversion from g    to Gg  (further x rho0)
 
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: btmsk   ! T-point basin interior masks
@@ -346,6 +346,6 @@
          IF( dia_ptr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ptr_init : unable to allocate arrays' )
 
-         rc_pwatt = rc_pwatt * rau0_rcp          ! conversion from K.s-1 to PetaWatt
-         rc_ggram = rc_ggram * rau0              ! conversion from m3/s to Gg/s
+         rc_pwatt = rc_pwatt * rho0_rcp          ! conversion from K.s-1 to PetaWatt
+         rc_ggram = rc_ggram * rho0              ! conversion from m3/s to Gg/s
 
          IF( lk_mpp )   CALL mpp_ini_znl( numout )     ! Define MPI communicator for zonal sum
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DIA/diawri.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DIA/diawri.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DIA/diawri.F90	(revision 12443)
@@ -173,5 +173,5 @@
 
       IF ( iom_use("taubot") ) THEN                ! bottom stress
-         zztmp = rau0 * 0.25
+         zztmp = rho0 * 0.25
          z2d(:,:) = 0._wp
          DO_2D_00_00
@@ -212,5 +212,5 @@
       IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value
          ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
-         z2d(:,:) = rau0 * e1e2t(:,:)
+         z2d(:,:) = rho0 * e1e2t(:,:)
          DO jk = 1, jpk
             z3d(:,:,jk) = ww(:,:,jk) * z2d(:,:)
@@ -249,5 +249,5 @@
             z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) * tmask(ji,jj,jk)
          END_3D
-         CALL iom_put( "heatc", rau0_rcp * z2d )   ! vertically integrated heat content (J/m2)
+         CALL iom_put( "heatc", rho0_rcp * z2d )   ! vertically integrated heat content (J/m2)
       ENDIF
 
@@ -257,5 +257,5 @@
             z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk)
          END_3D
-         CALL iom_put( "saltc", rau0 * z2d )          ! vertically integrated salt content (PSU*kg/m2)
+         CALL iom_put( "saltc", rho0 * z2d )          ! vertically integrated salt content (PSU*kg/m2)
       ENDIF
       !
@@ -279,5 +279,5 @@
          z2d(:,:) = 0.e0
          DO jk = 1, jpkm1
-            z3d(:,:,jk) = rau0 * uu(:,:,jk,Kmm) * e2u(:,:) * e3u(:,:,jk,Kmm) * umask(:,:,jk)
+            z3d(:,:,jk) = rho0 * uu(:,:,jk,Kmm) * e2u(:,:) * e3u(:,:,jk,Kmm) * umask(:,:,jk)
             z2d(:,:) = z2d(:,:) + z3d(:,:,jk)
          END DO
@@ -308,5 +308,5 @@
          z3d(:,:,jpk) = 0.e0
          DO jk = 1, jpkm1
-            z3d(:,:,jk) = rau0 * vv(:,:,jk,Kmm) * e1v(:,:) * e3v(:,:,jk,Kmm) * vmask(:,:,jk)
+            z3d(:,:,jk) = rho0 * vv(:,:,jk,Kmm) * e1v(:,:) * e3v(:,:,jk,Kmm) * vmask(:,:,jk)
          END DO
          CALL iom_put( "v_masstr", z3d )              ! mass transport in j-direction
@@ -337,5 +337,5 @@
          END_3D
          CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
-         CALL iom_put( "tosmint", rau0 * z2d )        ! Vertical integral of temperature
+         CALL iom_put( "tosmint", rho0 * z2d )        ! Vertical integral of temperature
       ENDIF
       IF( iom_use("somint") ) THEN
@@ -345,5 +345,5 @@
          END_3D
          CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
-         CALL iom_put( "somint", rau0 * z2d )         ! Vertical integral of salinity
+         CALL iom_put( "somint", rho0 * z2d )         ! Vertical integral of salinity
       ENDIF
 
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DOM/dom_oce.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DOM/dom_oce.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DOM/dom_oce.F90	(revision 12443)
@@ -33,5 +33,5 @@
    LOGICAL , PUBLIC ::   ln_linssh      !: =T  linear free surface ==>> model level are fixed in time
    LOGICAL , PUBLIC ::   ln_meshmask    !: =T  create a mesh-mask file (mesh_mask.nc)
-   REAL(wp), PUBLIC ::   rn_Dt         !: time step for the dynamics and tracer
+   REAL(wp), PUBLIC ::   rn_Dt          !: time step for the dynamics and tracer
    REAL(wp), PUBLIC ::   rn_atfp        !: asselin time filter parameter
    LOGICAL , PUBLIC ::   ln_1st_euler   !: =T start with forward time step or not (=F)
@@ -49,18 +49,15 @@
    LOGICAL,  PUBLIC :: ln_bt_auto       !: Set number of barotropic iterations automatically
    INTEGER,  PUBLIC :: nn_bt_flt        !: Filter choice
-   INTEGER,  PUBLIC :: nn_baro          !: Number of barotropic iterations during one baroclinic step (rn_Dt)
+   INTEGER,  PUBLIC :: nn_e          !: Number of barotropic iterations during one baroclinic step (rn_Dt)
    REAL(wp), PUBLIC :: rn_bt_cmax       !: Maximum allowed courant number (used if ln_bt_auto=T)
    REAL(wp), PUBLIC :: rn_bt_alpha      !: Time stepping diffusion parameter
 
-
-   !                                   !! old non-DOCTOR names still used in the model
-   REAL(wp), PUBLIC ::   atfp           !: asselin time filter parameter
 
    !                                   !!! associated variables
    LOGICAL , PUBLIC ::   l_1st_euler    !: Euler 1st time-step flag (=T if ln_restart=F or ln_1st_euler=T)
    REAL(wp), PUBLIC ::   rDt, r1_Dt     !: Current model timestep and reciprocal
-                                        !:  = 2 * rn_Dt if leapfrog and l_1st_euler = F
-                                        !:  =     rn_Dt if leapfrog and l_1st_euler = T
-                                        !:  =     rn_Dt if RK3
+                                        !: rDt = 2 * rn_Dt if leapfrog and l_1st_euler = F
+                                        !:     =     rn_Dt if leapfrog and l_1st_euler = T
+                                        !:     =     rn_Dt if RK3
 
    !!----------------------------------------------------------------------
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DOM/domain.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DOM/domain.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DOM/domain.F90	(revision 12443)
@@ -408,11 +408,9 @@
          WRITE(numout,*) '      linear free surface (=T)                ln_linssh   = ', ln_linssh
          WRITE(numout,*) '      create mesh/mask file                   ln_meshmask = ', ln_meshmask
-         WRITE(numout,*) '      ocean time step                         rn_Dt      = ', rn_Dt
+         WRITE(numout,*) '      ocean time step                         rn_Dt       = ', rn_Dt
          WRITE(numout,*) '      asselin time filter parameter           rn_atfp     = ', rn_atfp
          WRITE(numout,*) '      online coarsening of dynamical fields   ln_crs      = ', ln_crs
       ENDIF
       !
-      !          ! conversion DOCTOR names into model names (this should disappear soon)
-      atfp = rn_atfp
       !! Initialise current model timestep rDt = 2*rn_Dt if MLF or rDt = rn_Dt if RK3
       rDt  = 2._wp * rn_Dt
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DOM/domvvl.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DOM/domvvl.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DOM/domvvl.F90	(revision 12443)
@@ -621,5 +621,5 @@
          ELSE
             tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 
-            &         + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
+            &         + rn_atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
          ENDIF
          tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DOM/phycst.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DOM/phycst.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DOM/phycst.F90	(revision 12443)
@@ -39,10 +39,10 @@
    REAL(wp), PUBLIC ::   rt0      = 273.15_wp        !: freezing point of fresh water [Kelvin]
 
-   REAL(wp), PUBLIC ::   rau0                        !: volumic mass of reference     [kg/m3]
-   REAL(wp), PUBLIC ::   r1_rau0                     !: = 1. / rau0                   [m3/kg]
+   REAL(wp), PUBLIC ::   rho0                        !: volumic mass of reference     [kg/m3]
+   REAL(wp), PUBLIC ::   r1_rho0                     !: = 1. / rho0                   [m3/kg]
    REAL(wp), PUBLIC ::   rcp                         !: ocean specific heat           [J/Kelvin]
    REAL(wp), PUBLIC ::   r1_rcp                      !: = 1. / rcp                    [Kelvin/J]
-   REAL(wp), PUBLIC ::   rau0_rcp                    !: = rau0 * rcp 
-   REAL(wp), PUBLIC ::   r1_rau0_rcp                 !: = 1. / ( rau0 * rcp )
+   REAL(wp), PUBLIC ::   rho0_rcp                    !: = rho0 * rcp 
+   REAL(wp), PUBLIC ::   r1_rho0_rcp                 !: = 1. / ( rho0 * rcp )
 
    REAL(wp), PUBLIC ::   emic     =    0.97_wp       !: emissivity of snow or ice (not used?)
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DYN/dynatf.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DYN/dynatf.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DYN/dynatf.F90	(revision 12443)
@@ -87,5 +87,5 @@
       !!             arrays to start the next time step:
       !!                (puu(Kmm),pvv(Kmm)) = (puu(Kmm),pvv(Kmm)) 
-      !!                                    + atfp [ (puu(Kbb),pvv(Kbb)) + (puu(Kaa),pvv(Kaa)) - 2 (puu(Kmm),pvv(Kmm)) ]
+      !!                                    + rn_atfp [ (puu(Kbb),pvv(Kbb)) + (puu(Kaa),pvv(Kaa)) - 2 (puu(Kmm),pvv(Kmm)) ]
       !!             Note that with flux form advection and non linear free surface,
       !!             the time filter is applied on thickness weighted velocity.
@@ -181,6 +181,6 @@
             !                             ! =============!
             DO_3D_11_11( 1, jpkm1 )
-               puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Kmm) + atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Kaa) )
-               pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Kmm) + atfp * ( pvv(ji,jj,jk,Kbb) - 2._wp * pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk,Kaa) )
+               puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Kmm) + rn_atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Kaa) )
+               pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Kmm) + rn_atfp * ( pvv(ji,jj,jk,Kbb) - 2._wp * pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk,Kaa) )
             END_3D
             !                             ! ================!
@@ -191,9 +191,9 @@
             ALLOCATE( ze3t_f(jpi,jpj,jpk), zwfld(jpi,jpj) )
             DO jk = 1, jpkm1
-               ze3t_f(:,:,jk) = pe3t(:,:,jk,Kmm) + atfp * ( pe3t(:,:,jk,Kbb) - 2._wp * pe3t(:,:,jk,Kmm) + pe3t(:,:,jk,Kaa) )
+               ze3t_f(:,:,jk) = pe3t(:,:,jk,Kmm) + rn_atfp * ( pe3t(:,:,jk,Kbb) - 2._wp * pe3t(:,:,jk,Kmm) + pe3t(:,:,jk,Kaa) )
             END DO
             ! Add volume filter correction: compatibility with tracer advection scheme
             ! => time filter + conservation correction
-            zcoef = atfp * rn_Dt * r1_rau0
+            zcoef = rn_atfp * rn_Dt * r1_rho0
             zwfld(:,:) = emp_b(:,:) - emp(:,:)
             IF ( ln_rnf ) zwfld(:,:) =  zwfld(:,:) - ( rnf_b(:,:) - rnf(:,:) )
@@ -207,5 +207,5 @@
             !     to manage rnf, isf and possibly in the futur icb, tide water glacier (...)
             !     ...(kt, coef, ktop, kbot, hz, fwf_b, fwf)
-            IF ( ln_isf ) CALL isf_dynatf( kt, Kmm, ze3t_f, atfp * rn_Dt )
+            IF ( ln_isf ) CALL isf_dynatf( kt, Kmm, ze3t_f, rn_atfp * rn_Dt )
             !
             pe3t(:,:,1:jpkm1,Kmm) = ze3t_f(:,:,1:jpkm1)        ! filtered scale factor at T-points
@@ -216,6 +216,6 @@
                CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), pe3v(:,:,:,Kmm), 'V' )
                DO_3D_11_11( 1, jpkm1 )
-                  puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Kmm) + atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Kaa) )
-                  pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Kmm) + atfp * ( pvv(ji,jj,jk,Kbb) - 2._wp * pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk,Kaa) )
+                  puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Kmm) + rn_atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Kaa) )
+                  pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Kmm) + rn_atfp * ( pvv(ji,jj,jk,Kbb) - 2._wp * pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk,Kaa) )
                END_3D
                !
@@ -234,6 +234,6 @@
                   zve3b = pe3v(ji,jj,jk,Kbb) * pvv(ji,jj,jk,Kbb)
                   !
-                  puu(ji,jj,jk,Kmm) = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
-                  pvv(ji,jj,jk,Kmm) = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
+                  puu(ji,jj,jk,Kmm) = ( zue3n + rn_atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
+                  pvv(ji,jj,jk,Kmm) = ( zve3n + rn_atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
                END_3D
                pe3u(:,:,1:jpkm1,Kmm) = ze3u_f(:,:,1:jpkm1)  
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DYN/dynspg.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DYN/dynspg.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DYN/dynspg.F90	(revision 12443)
@@ -67,6 +67,6 @@
       !!              ln_apr_dyn=T : the atmospheric pressure forcing is applied 
       !!             as the gradient of the inverse barometer ssh:
-      !!                apgu = - 1/rau0 di[apr] = 0.5*grav di[ssh_ib+ssh_ibb]
-      !!                apgv = - 1/rau0 dj[apr] = 0.5*grav dj[ssh_ib+ssh_ibb]
+      !!                apgu = - 1/rho0 di[apr] = 0.5*grav di[ssh_ib+ssh_ibb]
+      !!                apgv = - 1/rho0 dj[apr] = 0.5*grav dj[ssh_ib+ssh_ibb]
       !!             Note that as all external forcing a time averaging over a two rn_Dt
       !!             period is used to prevent the divergence of odd and even time step.
@@ -78,5 +78,5 @@
       !
       INTEGER  ::   ji, jj, jk                   ! dummy loop indices
-      REAL(wp) ::   z2dt, zg_2, zintp, zgrau0r, zld   ! local scalars
+      REAL(wp) ::   z2dt, zg_2, zintp, zgrho0r, zld   ! local scalars
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zpice
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv
@@ -123,5 +123,5 @@
             !
             IF (ln_scal_load) THEN
-               zld = rn_scal_load * grav
+               zld = rn_load * grav
                DO_2D_00_00
                   spgu(ji,jj) = spgu(ji,jj) + zld * ( pssh(ji+1,jj,Kmm) - pssh(ji,jj,Kmm) ) * r1_e1u(ji,jj)
@@ -134,6 +134,6 @@
             ALLOCATE( zpice(jpi,jpj) )
             zintp = REAL( MOD( kt-1, nn_fsbc ) ) / REAL( nn_fsbc )
-            zgrau0r     = - grav * r1_rau0
-            zpice(:,:) = (  zintp * snwice_mass(:,:) + ( 1.- zintp ) * snwice_mass_b(:,:)  ) * zgrau0r
+            zgrho0r     = - grav * r1_rho0
+            zpice(:,:) = (  zintp * snwice_mass(:,:) + ( 1.- zintp ) * snwice_mass_b(:,:)  ) * zgrho0r
             DO_2D_00_00
                spgu(ji,jj) = spgu(ji,jj) + ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
@@ -183,5 +183,5 @@
       NAMELIST/namdyn_spg/ ln_dynspg_exp       , ln_dynspg_ts,   &
       &                    ln_bt_fw, ln_bt_av  , ln_bt_auto  ,   &
-      &                    nn_baro , rn_bt_cmax, nn_bt_flt, rn_bt_alpha
+      &                    nn_e , rn_bt_cmax, nn_bt_flt, rn_bt_alpha
       !!----------------------------------------------------------------------
       !
@@ -222,5 +222,5 @@
       !
       IF( nspg == np_TS ) THEN   ! split-explicit scheme initialisation
-         CALL dyn_spg_ts_init          ! do it first: set nn_baro used to allocate some arrays later on
+         CALL dyn_spg_ts_init          ! do it first: set nn_e used to allocate some arrays later on
       ENDIF
       !
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DYN/dynspg_exp.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DYN/dynspg_exp.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DYN/dynspg_exp.F90	(revision 12443)
@@ -49,6 +49,6 @@
       !!              momentum trend the surface pressure gradient :
       !!                      (uu(rhs),vv(rhs)) = (uu(rhs),vv(rhs)) + (spgu,spgv)
-      !!              where spgu = -1/rau0 d/dx(ps) = -g/e1u di( ssh(now) )
-      !!                    spgv = -1/rau0 d/dy(ps) = -g/e2v dj( ssh(now) )
+      !!              where spgu = -1/rho0 d/dx(ps) = -g/e1u di( ssh(now) )
+      !!                    spgv = -1/rho0 d/dy(ps) = -g/e2v dj( ssh(now) )
       !!
       !! ** Action :   (puu(:,:,:,Krhs),pvv(:,:,:,Krhs))   trend of horizontal velocity increased by 
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DYN/dynspg_ts.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DYN/dynspg_ts.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DYN/dynspg_ts.F90	(revision 12443)
@@ -72,6 +72,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   un_adv , vn_adv   !: Advection vel. at "now" barocl. step
    !
-   INTEGER, SAVE :: icycle      ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro
-   REAL(wp),SAVE :: rdtbt       ! Barotropic time step
+   INTEGER, SAVE :: icycle      ! Number of barotropic sub-steps for each internal step nn_e <= 2.5 nn_e
+   REAL(wp),SAVE :: rDt_e       ! Barotropic time step
    !
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)   ::   wgtbtp1, wgtbtp2   ! 1st & 2nd weights used in time filtering of barotropic fields
@@ -102,5 +102,5 @@
       ierr(:) = 0
       !
-      ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT=ierr(1) )
+      ALLOCATE( wgtbtp1(3*nn_e), wgtbtp2(3*nn_e), zwz(jpi,jpj), STAT=ierr(1) )
       IF( ln_dynvor_een .OR. ln_dynvor_eeT )   &
          &     ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , ftsw(jpi,jpj) , ftse(jpi,jpj), STAT=ierr(2)   )
@@ -185,5 +185,5 @@
       ll_fw_start = .FALSE.
       !                                         ! time offset in steps for bdy data update
-      IF( .NOT.ln_bt_fw ) THEN   ;   noffset = - nn_baro
+      IF( .NOT.ln_bt_fw ) THEN   ;   noffset = - nn_e
       ELSE                       ;   noffset =   0 
       ENDIF
@@ -303,9 +303,9 @@
       IF( ln_bt_fw ) THEN                        ! Add wind forcing
          DO_2D_00_00
-            zu_frc(ji,jj) =  zu_frc(ji,jj) + r1_rau0 * utau(ji,jj) * r1_hu(ji,jj,Kmm)
-            zv_frc(ji,jj) =  zv_frc(ji,jj) + r1_rau0 * vtau(ji,jj) * r1_hv(ji,jj,Kmm)
+            zu_frc(ji,jj) =  zu_frc(ji,jj) + r1_rho0 * utau(ji,jj) * r1_hu(ji,jj,Kmm)
+            zv_frc(ji,jj) =  zv_frc(ji,jj) + r1_rho0 * vtau(ji,jj) * r1_hv(ji,jj,Kmm)
          END_2D
       ELSE
-         zztmp = r1_rau0 * r1_2
+         zztmp = r1_rho0 * r1_2
          DO_2D_00_00
             zu_frc(ji,jj) =  zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu(ji,jj,Kmm)
@@ -320,7 +320,7 @@
       !                                   ! ---------------------------------------------------  !
       IF (ln_bt_fw) THEN                          ! FORWARD integration: use kt+1/2 fluxes (NOW+1/2)
-         zssh_frc(:,:) = r1_rau0 * ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) )
+         zssh_frc(:,:) = r1_rho0 * ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) )
       ELSE                                        ! CENTRED integration: use kt-1/2 + kt+1/2 fluxes (NOW)
-         zztmp = r1_rau0 * r1_2
+         zztmp = r1_rho0 * r1_2
          zssh_frc(:,:) = zztmp * (  emp(:,:)        + emp_b(:,:)                    &
                 &                 - rnf(:,:)        - rnf_b(:,:)                    &
@@ -429,5 +429,5 @@
          ! Update tide potential at the beginning of current time substep
          IF( ln_tide_pot .AND. ln_tide ) THEN
-            zt0substep = REAL(nsec_day, wp) - 0.5_wp*rn_Dt + (jn + noffset - 1) * rn_Dt / REAL(nn_baro, wp)
+            zt0substep = REAL(nsec_day, wp) - 0.5_wp*rn_Dt + (jn + noffset - 1) * rn_Dt / REAL(nn_e, wp)
             CALL upd_tide(zt0substep, Kmm)
          END IF
@@ -495,5 +495,5 @@
          IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) CALL agrif_dyn_ts_flux( jn, zhU, zhV )
 #endif
-         IF( ln_wd_il )   CALL wad_lmt_bt(zhU, zhV, sshn_e, zssh_frc, rdtbt)    !!gm wad_lmt_bt use of lbc_lnk on zhU, zhV
+         IF( ln_wd_il )   CALL wad_lmt_bt(zhU, zhV, sshn_e, zssh_frc, rDt_e)    !!gm wad_lmt_bt use of lbc_lnk on zhU, zhV
 
          IF( ln_wd_dl ) THEN           ! un_e and vn_e are set to zero at faces where 
@@ -510,5 +510,5 @@
          DO_2D_00_00
             zhdiv = (   zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1)   ) * r1_e1e2t(ji,jj)
-            ssha_e(ji,jj) = (  sshn_e(ji,jj) - rdtbt * ( zssh_frc(ji,jj) + zhdiv )  ) * ssmask(ji,jj)
+            ssha_e(ji,jj) = (  sshn_e(ji,jj) - rDt_e * ( zssh_frc(ji,jj) + zhdiv )  ) * ssmask(ji,jj)
          END_2D
          !
@@ -552,5 +552,5 @@
          !
          !                             ! Surface pressure gradient
-         zldg = ( 1._wp - rn_scal_load ) * grav    ! local factor
+         zldg = ( 1._wp - rn_load ) * grav    ! local factor
          DO_2D_00_00
             zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
@@ -600,5 +600,5 @@
             DO_2D_00_00
                ua_e(ji,jj) = (                                 un_e(ji,jj)   & 
-                         &     + rdtbt * (                   zu_spg(ji,jj)   &
+                         &     + rDt_e * (                   zu_spg(ji,jj)   &
                          &                                 + zu_trd(ji,jj)   &
                          &                                 + zu_frc(ji,jj) ) & 
@@ -606,5 +606,5 @@
 
                va_e(ji,jj) = (                                 vn_e(ji,jj)   &
-                         &     + rdtbt * (                   zv_spg(ji,jj)   &
+                         &     + rDt_e * (                   zv_spg(ji,jj)   &
                          &                                 + zv_trd(ji,jj)   &
                          &                                 + zv_frc(ji,jj) ) &
@@ -625,10 +625,10 @@
                !
                ua_e(ji,jj) = (               hu_e  (ji,jj) *   un_e (ji,jj)      & 
-                    &            + rdtbt * (  zhu_bck        * zu_spg (ji,jj)  &   !
+                    &            + rDt_e * (  zhu_bck        * zu_spg (ji,jj)  &   !
                     &                       + zhup2_e(ji,jj) * zu_trd (ji,jj)  &   !
                     &                       +  hu(ji,jj,Kmm) * zu_frc (ji,jj)  )   ) * z1_hu
                !
                va_e(ji,jj) = (               hv_e  (ji,jj) *   vn_e (ji,jj)      &
-                    &            + rdtbt * (  zhv_bck        * zv_spg (ji,jj)  &   !
+                    &            + rDt_e * (  zhv_bck        * zv_spg (ji,jj)  &   !
                     &                       + zhvp2_e(ji,jj) * zv_trd (ji,jj)  &   !
                     &                       +  hv(ji,jj,Kmm) * zv_frc (ji,jj)  )   ) * z1_hv
@@ -638,6 +638,6 @@
          IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs
             DO_2D_00_00
-                  ua_e(ji,jj) =  ua_e(ji,jj) /(1.0 -   rdtbt * zCdU_u(ji,jj) * hur_e(ji,jj))
-                  va_e(ji,jj) =  va_e(ji,jj) /(1.0 -   rdtbt * zCdU_v(ji,jj) * hvr_e(ji,jj))
+                  ua_e(ji,jj) =  ua_e(ji,jj) /(1.0 -   rDt_e * zCdU_u(ji,jj) * hur_e(ji,jj))
+                  va_e(ji,jj) =  va_e(ji,jj) /(1.0 -   rDt_e * zCdU_v(ji,jj) * hvr_e(ji,jj))
             END_2D
          ENDIF
@@ -707,9 +707,9 @@
                zvn_save = vn_adv(ji,jj)
                !                          ! apply the previously computed correction 
-               un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - atfp * un_bf(ji,jj) )
-               vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - atfp * vn_bf(ji,jj) )
+               un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - rn_atfp * un_bf(ji,jj) )
+               vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - rn_atfp * vn_bf(ji,jj) )
                !                          ! Update corrective fluxes for next time step
-               un_bf(ji,jj)  = atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) )
-               vn_bf(ji,jj)  = atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) )
+               un_bf(ji,jj)  = rn_atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) )
+               vn_bf(ji,jj)  = rn_atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) )
                !                          ! Save integrated transport for next computation
                ub2_b(ji,jj) = zun_save
@@ -809,5 +809,5 @@
       LOGICAL, INTENT(in) ::   ll_fw      ! forward time splitting =.true.
       INTEGER, INTENT(inout) :: jpit      ! cycle length    
-      REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) ::   zwgt1, & ! Primary weights
+      REAL(wp), DIMENSION(3*nn_e), INTENT(inout) ::   zwgt1, & ! Primary weights
                                                          zwgt2    ! Secondary weights
       
@@ -821,7 +821,7 @@
       ! Set time index when averaged value is requested
       IF (ll_fw) THEN 
-         jic = nn_baro
+         jic = nn_e
       ELSE
-         jic = 2 * nn_baro
+         jic = 2 * nn_e
       ENDIF
 
@@ -829,5 +829,5 @@
       IF (ll_av) THEN
            ! Define simple boxcar window for primary weights 
-           ! (width = nn_baro, centered around jic)     
+           ! (width = nn_e, centered around jic)     
          SELECT CASE ( nn_bt_flt )
               CASE( 0 )  ! No averaging
@@ -835,7 +835,7 @@
                  jpit = jic
 
-              CASE( 1 )  ! Boxcar, width = nn_baro
-                 DO jn = 1, 3*nn_baro
-                    za1 = ABS(float(jn-jic))/float(nn_baro) 
+              CASE( 1 )  ! Boxcar, width = nn_e
+                 DO jn = 1, 3*nn_e
+                    za1 = ABS(float(jn-jic))/float(nn_e) 
                     IF (za1 < 0.5_wp) THEN
                       zwgt1(jn) = 1._wp
@@ -844,7 +844,7 @@
                  ENDDO
 
-              CASE( 2 )  ! Boxcar, width = 2 * nn_baro
-                 DO jn = 1, 3*nn_baro
-                    za1 = ABS(float(jn-jic))/float(nn_baro) 
+              CASE( 2 )  ! Boxcar, width = 2 * nn_e
+                 DO jn = 1, 3*nn_e
+                    za1 = ABS(float(jn-jic))/float(nn_e) 
                     IF (za1 < 1._wp) THEN
                       zwgt1(jn) = 1._wp
@@ -976,8 +976,8 @@
 
       ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
-      IF( ln_bt_auto )   nn_baro = CEILING( rn_Dt / rn_bt_cmax * zcmax)
+      IF( ln_bt_auto )   nn_e = CEILING( rn_Dt / rn_bt_cmax * zcmax)
       
-      rdtbt = rn_Dt / REAL( nn_baro , wp )
-      zcmax = zcmax * rdtbt
+      rDt_e = rn_Dt / REAL( nn_e , wp )
+      zcmax = zcmax * rDt_e
       ! Print results
       IF(lwp) WRITE(numout,*)
@@ -985,12 +985,12 @@
       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
       IF( ln_bt_auto ) THEN
-         IF(lwp) WRITE(numout,*) '     ln_ts_auto =.true. Automatically set nn_baro '
+         IF(lwp) WRITE(numout,*) '     ln_ts_auto =.true. Automatically set nn_e '
          IF(lwp) WRITE(numout,*) '     Max. courant number allowed: ', rn_bt_cmax
       ELSE
-         IF(lwp) WRITE(numout,*) '     ln_ts_auto=.false.: Use nn_baro in namelist   nn_baro = ', nn_baro
+         IF(lwp) WRITE(numout,*) '     ln_ts_auto=.false.: Use nn_e in namelist   nn_e = ', nn_e
       ENDIF
 
       IF(ln_bt_av) THEN
-         IF(lwp) WRITE(numout,*) '     ln_bt_av =.true.  ==> Time averaging over nn_baro time steps is on '
+         IF(lwp) WRITE(numout,*) '     ln_bt_av =.true.  ==> Time averaging over nn_e time steps is on '
       ELSE
          IF(lwp) WRITE(numout,*) '     ln_bt_av =.false. => No time averaging of barotropic variables '
@@ -1012,12 +1012,12 @@
       SELECT CASE ( nn_bt_flt )
          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '           Dirac'
-         CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = nn_baro'
-         CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = 2*nn_baro' 
+         CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = nn_e'
+         CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = 2*nn_e' 
          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' )
       END SELECT
       !
       IF(lwp) WRITE(numout,*) ' '
-      IF(lwp) WRITE(numout,*) '     nn_baro = ', nn_baro
-      IF(lwp) WRITE(numout,*) '     Barotropic time step [s] is :', rdtbt
+      IF(lwp) WRITE(numout,*) '     nn_e = ', nn_e
+      IF(lwp) WRITE(numout,*) '     Barotropic time step [s] is :', rDt_e
       IF(lwp) WRITE(numout,*) '     Maximum Courant number is   :', zcmax
       !
@@ -1031,5 +1031,5 @@
       ENDIF
       IF( zcmax>0.9_wp ) THEN
-         CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )          
+         CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_e !' )          
       ENDIF
       !
@@ -1430,5 +1430,5 @@
       !
       IF( ln_wd_il ) THEN      ! W/D : use the "clipped" bottom friction   !!gm   explain WHY, please !
-         zztmp = -1._wp / rdtbt
+         zztmp = -1._wp / rDt_e
          DO_2D_00_00
             pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) *  wdrampu(ji,jj) * MAX(                                 & 
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DYN/dynzdf.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DYN/dynzdf.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DYN/dynzdf.F90	(revision 12443)
@@ -261,5 +261,5 @@
          ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) + r_vvl * e3u(ji,jj,1,Kaa) 
          puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + rDt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
-            &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1) 
+            &                                      / ( ze3ua * rho0 ) * umask(ji,jj,1) 
       END_2D
       DO_3D_00_00( 2, jpkm1 )
@@ -385,5 +385,5 @@
          ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) + r_vvl * e3v(ji,jj,1,Kaa) 
          pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + rDt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
-            &                                      / ( ze3va * rau0 ) * vmask(ji,jj,1) 
+            &                                      / ( ze3va * rho0 ) * vmask(ji,jj,1) 
       END_2D
       DO_3D_00_00( 2, jpkm1 )
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DYN/sshwzv.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DYN/sshwzv.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DYN/sshwzv.F90	(revision 12443)
@@ -88,5 +88,5 @@
       ENDIF
       !
-      zcoef = 0.5_wp * r1_rau0
+      zcoef = 0.5_wp * r1_rho0
 
       !                                           !------------------------------!
@@ -222,6 +222,6 @@
       !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
       !!              from the filter, see Leclair and Madec 2010) and swap :
-      !!                pssh(:,:,Kmm) = pssh(:,:,Kaa) + atfp * ( pssh(:,:,Kbb) -2 pssh(:,:,Kmm) + pssh(:,:,Kaa) )
-      !!                            - atfp * rn_Dt * ( emp_b - emp ) / rau0
+      !!                pssh(:,:,Kmm) = pssh(:,:,Kaa) + rn_atfp * ( pssh(:,:,Kbb) -2 pssh(:,:,Kmm) + pssh(:,:,Kaa) )
+      !!                            - rn_atfp * rn_Dt * ( emp_b - emp ) / rho0
       !!
       !! ** action  : - pssh(:,:,Kmm) time filtered
@@ -246,7 +246,7 @@
       IF ( .NOT.( l_1st_euler ) ) THEN   ! Only do time filtering for leapfrog timesteps
          !                                                  ! filtered "now" field
-         pssh(:,:,Kmm) = pssh(:,:,Kmm) + atfp * ( pssh(:,:,Kbb) - 2 * pssh(:,:,Kmm) + pssh(:,:,Kaa) )
+         pssh(:,:,Kmm) = pssh(:,:,Kmm) + rn_atfp * ( pssh(:,:,Kbb) - 2 * pssh(:,:,Kmm) + pssh(:,:,Kaa) )
          IF( .NOT.ln_linssh ) THEN                          ! "now" <-- with forcing removed
-            zcoef = atfp * rn_Dt * r1_rau0
+            zcoef = rn_atfp * rn_Dt * r1_rho0
             pssh(:,:,Kmm) = pssh(:,:,Kmm) - zcoef * (     emp_b(:,:) - emp   (:,:)   &
                &                             - rnf_b(:,:)        + rnf   (:,:)       &
@@ -255,5 +255,5 @@
 
             ! ice sheet coupling
-            IF ( ln_isf .AND. ln_isfcpl .AND. kt == nit000+1) pssh(:,:,Kbb) = pssh(:,:,Kbb) - atfp * rn_Dt * ( risfcpl_ssh(:,:) - 0.0 ) * ssmask(:,:)
+            IF ( ln_isf .AND. ln_isfcpl .AND. kt == nit000+1) pssh(:,:,Kbb) = pssh(:,:,Kbb) - rn_atfp * rn_Dt * ( risfcpl_ssh(:,:) - 0.0 ) * ssmask(:,:)
 
          ENDIF
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DYN/wet_dry.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DYN/wet_dry.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DYN/wet_dry.F90	(revision 12443)
@@ -270,5 +270,5 @@
 
 
-   SUBROUTINE wad_lmt_bt( zflxu, zflxv, sshn_e, zssh_frc, rdtbt )
+   SUBROUTINE wad_lmt_bt( zflxu, zflxv, sshn_e, zssh_frc, rDt_e )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE wad_lmt  ***
@@ -280,5 +280,5 @@
       !! ** Action  : - calculate flux limiter and W/D flag
       !!----------------------------------------------------------------------
-      REAL(wp)                , INTENT(in   ) ::   rdtbt    ! ocean time-step index
+      REAL(wp)                , INTENT(in   ) ::   rDt_e    ! ocean time-step index
       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zflxu,  zflxv, sshn_e, zssh_frc  
       !
@@ -299,5 +299,5 @@
       zdepwd = 50._wp   ! maximum depth that ocean cells can have W/D processes
       !
-      z2dt = rdtbt   
+      z2dt = rDt_e   
       !
       zflxp(:,:)   = 0._wp
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/FLO/flowri.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/FLO/flowri.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/FLO/flowri.F90	(revision 12443)
@@ -122,5 +122,5 @@
                ztem(jfl) = ts(iafloc,ibfloc,icfl,jp_tem,Kmm)
                zsal (jfl) = ts(iafloc,ibfloc,icfl,jp_sal,Kmm)
-               zrho (jfl) = (rhd(iafloc,ibfloc,icfl)+1)*rau0
+               zrho (jfl) = (rhd(iafloc,ibfloc,icfl)+1)*rho0
 
             ENDIF
@@ -142,5 +142,5 @@
             ztem(jfl) = ts(iafloc,ibfloc,icfl,jp_tem,Kmm)
             zsal(jfl) = ts(iafloc,ibfloc,icfl,jp_sal,Kmm)
-            zrho(jfl) = (rhd(iafloc,ibfloc,icfl)+1)*rau0
+            zrho(jfl) = (rhd(iafloc,ibfloc,icfl)+1)*rho0
           
          ENDIF
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/IOM/restart.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/IOM/restart.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/IOM/restart.F90	(revision 12443)
@@ -262,6 +262,6 @@
       IF ( ln_diurnal_only ) THEN 
          IF(lwp) WRITE( numout, * ) &
-         &   "rst_read:- ln_diurnal_only set, setting rhop=rau0" 
-         rhop = rau0
+         &   "rst_read:- ln_diurnal_only set, setting rhop=rho0" 
+         rhop = rho0
          CALL iom_get( numror, jpdom_autoglo, 'tn'     , w3d, ldxios = lrxios ) 
          ts(:,:,1,jp_tem,Kmm) = w3d(:,:,1)
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ISF/isfcav.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ISF/isfcav.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ISF/isfcav.F90	(revision 12443)
@@ -24,5 +24,5 @@
    USE oce      , ONLY: ts                              ! ocean tracers
    USE par_oce  , ONLY: jpi,jpj                         ! ocean space and time domain
-   USE phycst   , ONLY: grav,rau0,rau0_rcp,r1_rau0_rcp  ! physical constants
+   USE phycst   , ONLY: grav,rho0,rho0_rcp,r1_rho0_rcp  ! physical constants
    USE eosbn2   , ONLY: ln_teos10                       ! use ln_teos10 or not
    !
@@ -85,5 +85,5 @@
       !
       ! initialisation
-      IF (TRIM(cn_gammablk) == 'vel_stab' ) zqoce_b (:,:) = ptsc(:,:,jp_tem) * rau0_rcp ! last time step total heat fluxes (to speed up convergence)
+      IF (TRIM(cn_gammablk) == 'vel_stab' ) zqoce_b (:,:) = ptsc(:,:,jp_tem) * rho0_rcp ! last time step total heat fluxes (to speed up convergence)
       !
       ! compute ice shelf melting
@@ -142,5 +142,5 @@
       !
       ! set temperature content
-      ptsc(:,:,jp_tem) = - zqh(:,:) * r1_rau0_rcp
+      ptsc(:,:,jp_tem) = - zqh(:,:) * r1_rho0_rcp
       !
       ! write restart variables (qoceisf, qhcisf, fwfisf for now and before)
@@ -215,9 +215,9 @@
             risf_lamb1 =-0.0564_wp
             risf_lamb2 = 0.0773_wp
-            risf_lamb3 =-7.8633e-8 * grav * rau0
+            risf_lamb3 =-7.8633e-8 * grav * rho0
          ELSE                  ! linearisation from table 4 (Asay-Davis et al., 2015)
             risf_lamb1 =-0.0573_wp
             risf_lamb2 = 0.0832_wp
-            risf_lamb3 =-7.5300e-8 * grav * rau0
+            risf_lamb3 =-7.5300e-8 * grav * rho0
          ENDIF
 
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ISF/isfcavmlt.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ISF/isfcavmlt.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ISF/isfcavmlt.F90	(revision 12443)
@@ -17,5 +17,5 @@
 
    USE dom_oce                            ! ocean space and time domain
-   USE phycst , ONLY: rcp, rau0, rau0_rcp ! physical constants
+   USE phycst , ONLY: rcp, rho0, rho0_rcp ! physical constants
    USE eosbn2 , ONLY: eos_fzp             ! equation of state
 
@@ -161,5 +161,5 @@
       !
       ! compute ocean-ice heat flux and then derive fwf assuming that ocean heat flux equal latent heat
-      pqfwf(:,:) = - pgt(:,:) * rau0_rcp * zthd(:,:) / rLfusisf    ! fresh water flux  ( > 0 out )
+      pqfwf(:,:) = - pgt(:,:) * rho0_rcp * zthd(:,:) / rLfusisf    ! fresh water flux  ( > 0 out )
       pqoce(:,:) = - pqfwf(:,:) * rLfusisf                         ! ocea-ice flux     ( > 0 out )
       pqhc (:,:) =   pqfwf(:,:) * ztfrz(:,:) * rcp                 ! heat content flux ( > 0 out )
@@ -213,6 +213,6 @@
          !
          ! compute coeficient to solve the 2nd order equation
-         zeps1 = rau0_rcp * pgt(ji,jj)
-         zeps2 = rLfusisf * rau0 * pgs(ji,jj)
+         zeps1 = rho0_rcp * pgt(ji,jj)
+         zeps2 = rLfusisf * rho0 * pgs(ji,jj)
          zeps3 = rhoisf * rcpisf * rkappa / MAX(risfdep(ji,jj),zeps)
          zeps4 = risf_lamb2 + risf_lamb3 * risfdep(ji,jj)
@@ -238,6 +238,6 @@
          !
          ! compute the upward water and heat flux (eq. 24 and eq. 26)
-         pqfwf(ji,jj) = rau0     * pgs(ji,jj) * ( zsfrz - pstbl(ji,jj) ) / MAX(zsfrz,zeps) ! fresh water flux    (> 0 out)
-         pqoce(ji,jj) = rau0_rcp * pgt(ji,jj) * zthd (ji,jj)                               ! ocean-ice heat flux (> 0 out)
+         pqfwf(ji,jj) = rho0     * pgs(ji,jj) * ( zsfrz - pstbl(ji,jj) ) / MAX(zsfrz,zeps) ! fresh water flux    (> 0 out)
+         pqoce(ji,jj) = rho0_rcp * pgt(ji,jj) * zthd (ji,jj)                               ! ocean-ice heat flux (> 0 out)
          pqhc (ji,jj) = rcp      * pqfwf(ji,jj) * ztfrz(ji,jj)                             ! heat content   flux (> 0 out)
          !
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ISF/isfdynatf.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ISF/isfdynatf.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ISF/isfdynatf.F90	(revision 12443)
@@ -13,5 +13,5 @@
    USE isf_oce
 
-   USE phycst , ONLY: r1_rau0         ! physical constant
+   USE phycst , ONLY: r1_rho0         ! physical constant
    USE dom_oce, ONLY: tmask, ssmask, ht, e3t, r1_e1e2t   ! time and space domain
 
@@ -39,5 +39,5 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pe3t_f   ! time filtered scale factor to be corrected
       !
-      REAL(wp)                        , INTENT(in   ) :: pcoef    ! atfp * rn_Dt * r1_rau0
+      REAL(wp)                        , INTENT(in   ) :: pcoef    ! rn_atfp * rn_Dt * r1_rho0
       !!--------------------------------------------------------------------
       INTEGER :: jk  ! loop index
@@ -70,5 +70,5 @@
       REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) :: pfrac, phtbl    ! fraction of bottom cell included in tbl, tbl thickness
       REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) :: pfwf , pfwf_b   ! now/before fwf
-      REAL(wp),                         INTENT(in   ) :: pcoef           ! atfp * rn_Dt * r1_rau0
+      REAL(wp),                         INTENT(in   ) :: pcoef           ! rn_atfp * rn_Dt * r1_rho0
       !!----------------------------------------------------------------------
       INTEGER :: ji,jj,jk
@@ -77,5 +77,5 @@
       !
       ! compute fwf conservation correction
-      zfwfinc(:,:) = pcoef * ( pfwf_b(:,:) - pfwf(:,:) ) / ( ht(:,:) + 1._wp - ssmask(:,:) ) * r1_rau0
+      zfwfinc(:,:) = pcoef * ( pfwf_b(:,:) - pfwf(:,:) ) / ( ht(:,:) + 1._wp - ssmask(:,:) ) * r1_rho0
       !
       ! add the increment
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ISF/isfhdiv.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ISF/isfhdiv.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ISF/isfhdiv.F90	(revision 12443)
@@ -16,5 +16,5 @@
 
    USE dom_oce                ! time and space domain
-   USE phycst , ONLY: r1_rau0 ! physical constant
+   USE phycst , ONLY: r1_rho0 ! physical constant
    USE in_out_manager         !
 
@@ -96,5 +96,5 @@
       !
       ! compute integrated divergence correction
-      zhdiv(:,:) = 0.5_wp * ( pfwf(:,:) + pfwf_b(:,:) ) * r1_rau0 / phtbl(:,:)
+      zhdiv(:,:) = 0.5_wp * ( pfwf(:,:) + pfwf_b(:,:) ) * r1_rho0 / phtbl(:,:)
       !
       ! update divergence at each level affected by ice shelf top boundary layer
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ISF/isfpar.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ISF/isfpar.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ISF/isfpar.F90	(revision 12443)
@@ -24,5 +24,5 @@
    USE dom_oce  , ONLY: bathy          ! ocean space and time domain
    USE par_oce  , ONLY: jpi,jpj        ! ocean space and time domain
-   USE phycst   , ONLY: r1_rau0_rcp    ! physical constants
+   USE phycst   , ONLY: r1_rho0_rcp    ! physical constants
    !
    USE in_out_manager ! I/O manager
@@ -88,5 +88,5 @@
       !
       ! set temperature content
-      ptsc(:,:,jp_tem) = zqh(:,:) * r1_rau0_rcp
+      ptsc(:,:,jp_tem) = zqh(:,:) * r1_rho0_rcp
       !
       ! write restart variables (qoceisf, qhcisf, fwfisf for now and before)
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ISF/isfparmlt.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ISF/isfparmlt.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ISF/isfparmlt.F90	(revision 12443)
@@ -13,5 +13,5 @@
    USE dom_oce                  ! ocean space and time domain
    USE oce    , ONLY: ts        ! ocean dynamics and tracers
-   USE phycst , ONLY: rcp, rau0 ! physical constants
+   USE phycst , ONLY: rcp, rho0 ! physical constants
    USE eosbn2 , ONLY: eos_fzp   ! equation of state
 
@@ -148,5 +148,5 @@
       !
       ! 2. ------------Net heat flux and fresh water flux due to the ice shelf
-      pqoce(:,:) =   rau0 * rcp * rn_gammat0 * risfLeff(:,:) * e1t(:,:) * ( ztavg(:,:) - ztfrz(:,:) ) * r1_e1e2t(:,:)
+      pqoce(:,:) =   rho0 * rcp * rn_gammat0 * risfLeff(:,:) * e1t(:,:) * ( ztavg(:,:) - ztfrz(:,:) ) * r1_e1e2t(:,:)
       pqfwf(:,:) = - pqoce(:,:) / rLfusisf             ! derived from the latent heat flux
       pqhc (:,:) =   pqfwf(:,:) * ztfrz(:,:) * rcp     ! heat content flux 
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/LDF/ldftra.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/LDF/ldftra.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/LDF/ldftra.F90	(revision 12443)
@@ -820,5 +820,5 @@
       !
       IF( iom_use('weiv_masstr') ) THEN   ! vertical mass transport & its square value
-         zw2d(:,:) = rau0 * e1e2t(:,:)
+         zw2d(:,:) = rho0 * e1e2t(:,:)
          DO jk = 1, jpk
             zw3d(:,:,jk) = zw3d(:,:,jk) * zw2d(:,:)
@@ -830,10 +830,10 @@
          zw3d(:,:,:) = 0.e0
          DO jk = 1, jpkm1
-            zw3d(:,:,jk) = rau0 * ( psi_uw(:,:,jk+1) - psi_uw(:,:,jk) ) 
+            zw3d(:,:,jk) = rho0 * ( psi_uw(:,:,jk+1) - psi_uw(:,:,jk) ) 
          END DO
          CALL iom_put( "ueiv_masstr", zw3d )                  ! mass transport in i-direction
       ENDIF
       !
-      zztmp = 0.5_wp * rau0 * rcp 
+      zztmp = 0.5_wp * rho0 * rcp 
       IF( iom_use('ueiv_heattr') .OR. iom_use('ueiv_heattr3d') ) THEN
         zw2d(:,:)   = 0._wp 
@@ -853,5 +853,5 @@
          zw3d(:,:,:) = 0.e0
          DO jk = 1, jpkm1
-            zw3d(:,:,jk) = rau0 * ( psi_vw(:,:,jk+1) - psi_vw(:,:,jk) ) 
+            zw3d(:,:,jk) = rho0 * ( psi_vw(:,:,jk+1) - psi_vw(:,:,jk) ) 
          END DO
          CALL iom_put( "veiv_masstr", zw3d )                  ! mass transport in i-direction
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/SBC/fldread.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/SBC/fldread.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/SBC/fldread.F90	(revision 12443)
@@ -172,5 +172,5 @@
       ! Note that all varibles starting by nsec_* are shifted time by +1/2 time step to be centrered
       IF( PRESENT(kit) ) THEN   ! ignore kn_fsbc in this case
-         isecsbc = nsec_year + nsec1jan000 + NINT( (     REAL(      kit,wp) + zt_offset ) * rn_Dt / REAL(nn_baro,wp) )
+         isecsbc = nsec_year + nsec1jan000 + NINT( (     REAL(      kit,wp) + zt_offset ) * rn_Dt / REAL(nn_e,wp) )
       ELSE                      ! middle of sbc time step
          ! note: we use kn_fsbc-1 because nsec_year is defined at the middle of the current time step
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/SBC/sbcapr.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/SBC/sbcapr.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/SBC/sbcapr.F90	(revision 12443)
@@ -36,5 +36,5 @@
    
    REAL(wp) ::   tarea                ! whole domain mean masked ocean surface
-   REAL(wp) ::   r1_grau              ! = 1.e0 / (grav * rau0)
+   REAL(wp) ::   r1_grau              ! = 1.e0 / (grav * rho0)
    
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_apr   ! structure of input fields (file informations, fields read)
@@ -98,5 +98,5 @@
       ENDIF
       !
-      r1_grau = 1.e0 / (grav * rau0)               !* constant for optimization
+      r1_grau = 1.e0 / (grav * rho0)               !* constant for optimization
       !
       !                                            !* control check
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/SBC/sbccpl.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/SBC/sbccpl.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/SBC/sbccpl.F90	(revision 12443)
@@ -193,5 +193,5 @@
 
    REAL(wp) ::   rpref = 101000._wp   ! reference atmospheric pressure[N/m2] 
-   REAL(wp) ::   r1_grau              ! = 1.e0 / (grav * rau0) 
+   REAL(wp) ::   r1_grau              ! = 1.e0 / (grav * rho0) 
 
    INTEGER , ALLOCATABLE, SAVE, DIMENSION(:) ::   nrcvinfo           ! OASIS info argument
@@ -1250,5 +1250,5 @@
           IF( kt /= nit000 )   ssh_ibb(:,:) = ssh_ib(:,:)    !* Swap of ssh_ib fields 
 
-          r1_grau = 1.e0 / (grav * rau0)               !* constant for optimization 
+          r1_grau = 1.e0 / (grav * rho0)               !* constant for optimization 
           ssh_ib(:,:) = - ( frcv(jpr_mslp)%z3(:,:,1) - rpref ) * r1_grau    ! equivalent ssh (inverse barometer) 
           apr   (:,:) =     frcv(jpr_mslp)%z3(:,:,1)                         !atmospheric pressure 
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/SBC/sbcfwb.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/SBC/sbcfwb.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/SBC/sbcfwb.F90	(revision 12443)
@@ -130,5 +130,5 @@
             a_fwb_b = a_fwb                           ! mean sea level taking into account the ice+snow
                                                       ! sum over the global domain
-            a_fwb   = glob_sum( 'sbcfwb', e1e2t(:,:) * ( ssh(:,:,Kmm) + snwice_mass(:,:) * r1_rau0 ) )
+            a_fwb   = glob_sum( 'sbcfwb', e1e2t(:,:) * ( ssh(:,:,Kmm) + snwice_mass(:,:) * r1_rho0 ) )
             a_fwb   = a_fwb * 1.e+3 / ( area * rday * 365. )     ! convert in Kg/m3/s = mm/s
 !!gm        !                                                      !!bug 365d year 
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/SBC/sbcice_cice.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/SBC/sbcice_cice.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/SBC/sbcice_cice.F90	(revision 12443)
@@ -13,5 +13,5 @@
    USE dom_oce         ! ocean space and time domain
    USE domvvl
-   USE phycst, only : rcp, rau0, r1_rau0, rhos, rhoi
+   USE phycst, only : rcp, rho0, r1_rho0, rhos, rhoi
    USE in_out_manager  ! I/O manager
    USE iom, ONLY : iom_put,iom_use              ! I/O manager library !!Joakim edit
@@ -228,6 +228,6 @@
       IF( .NOT.ln_rstart ) THEN
          IF( ln_ice_embd ) THEN            ! embedded sea-ice: deplete the initial ssh below sea-ice area
-            ssh(:,:,Kmm) = ssh(:,:,Kmm) - snwice_mass(:,:) * r1_rau0
-            ssh(:,:,Kbb) = ssh(:,:,Kbb) - snwice_mass(:,:) * r1_rau0
+            ssh(:,:,Kmm) = ssh(:,:,Kmm) - snwice_mass(:,:) * r1_rho0
+            ssh(:,:,Kbb) = ssh(:,:,Kbb) - snwice_mass(:,:) * r1_rho0
 
 !!gm This should be put elsewhere....   (same remark for limsbc)
@@ -417,5 +417,5 @@
 ! Freezing/melting potential
 ! Calculated over NEMO leapfrog timestep (hence 2*dt)
-      nfrzmlt(:,:) = rau0 * rcp * e3t_m(:,:) * ( Tocnfrz-sst_m(:,:) ) / ( 2.0*dt )
+      nfrzmlt(:,:) = rho0 * rcp * e3t_m(:,:) * ( Tocnfrz-sst_m(:,:) ) / ( 2.0*dt )
 
       ztmp(:,:) = nfrzmlt(:,:)
@@ -450,5 +450,5 @@
          zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
           !
-         zpice(:,:) = ssh_m(:,:) + (  zintn * snwice_mass(:,:) +  zintb * snwice_mass_b(:,:)  ) * r1_rau0
+         zpice(:,:) = ssh_m(:,:) + (  zintn * snwice_mass(:,:) +  zintb * snwice_mass_b(:,:)  ) * r1_rho0
           !
          !
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/SBC/sbcrnf.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/SBC/sbcrnf.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/SBC/sbcrnf.F90	(revision 12443)
@@ -137,18 +137,18 @@
          !                                                           ! set temperature & salinity content of runoffs
          IF( ln_rnf_tem ) THEN                                       ! use runoffs temperature data
-            rnf_tsc(:,:,jp_tem) = ( sf_t_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rau0
+            rnf_tsc(:,:,jp_tem) = ( sf_t_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rho0
             CALL eos_fzp( sss_m(:,:), ztfrz(:,:) )
             WHERE( sf_t_rnf(1)%fnow(:,:,1) == -999._wp )             ! if missing data value use SST as runoffs temperature
-               rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rau0
+               rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rho0
             END WHERE
          ELSE                                                        ! use SST as runoffs temperature
             !CEOD River is fresh water so must at least be 0 unless we consider ice
-            rnf_tsc(:,:,jp_tem) = MAX( sst_m(:,:), 0.0_wp ) * rnf(:,:) * r1_rau0
+            rnf_tsc(:,:,jp_tem) = MAX( sst_m(:,:), 0.0_wp ) * rnf(:,:) * r1_rho0
          ENDIF
          !                                                           ! use runoffs salinity data
-         IF( ln_rnf_sal )   rnf_tsc(:,:,jp_sal) = ( sf_s_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rau0
+         IF( ln_rnf_sal )   rnf_tsc(:,:,jp_sal) = ( sf_s_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rho0
          !                                                           ! else use S=0 for runoffs (done one for all in the init)
                                          CALL iom_put( 'runoffs'     , rnf(:,:)                         )   ! output runoff mass flux
-         IF( iom_use('hflx_rnf_cea') )   CALL iom_put( 'hflx_rnf_cea', rnf_tsc(:,:,jp_tem) * rau0 * rcp )   ! output runoff sensible heat (W/m2)
+         IF( iom_use('hflx_rnf_cea') )   CALL iom_put( 'hflx_rnf_cea', rnf_tsc(:,:,jp_tem) * rho0 * rcp )   ! output runoff sensible heat (W/m2)
       ENDIF
       !
@@ -210,5 +210,5 @@
             DO_2D_11_11
                DO jk = 1, nk_rnf(ji,jj)
-                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact * r1_rau0 / h_rnf(ji,jj)
+                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact * r1_rho0 / h_rnf(ji,jj)
                END DO
             END_2D
@@ -221,5 +221,5 @@
                !                          ! apply the runoff input flow
                DO jk = 1, nk_rnf(ji,jj)
-                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact * r1_rau0 / h_rnf(ji,jj)
+                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact * r1_rho0 / h_rnf(ji,jj)
                END DO
             END_2D
@@ -227,5 +227,5 @@
       ELSE                       !==   runoff put only at the surface   ==!
          h_rnf (:,:)   = e3t (:,:,1,Kmm)        ! update h_rnf to be depth of top box
-         phdivn(:,:,1) = phdivn(:,:,1) - ( rnf(:,:) + rnf_b(:,:) ) * zfact * r1_rau0 / e3t(:,:,1,Kmm)
+         phdivn(:,:,1) = phdivn(:,:,1) - ( rnf(:,:) + rnf_b(:,:) ) * zfact * r1_rho0 / e3t(:,:,1,Kmm)
       ENDIF
       !
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TDE/tide_mod.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TDE/tide_mod.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TDE/tide_mod.F90	(revision 12443)
@@ -64,5 +64,5 @@
    INTEGER , PUBLIC ::   nb_harmo            !: Number of active tidal components
    REAL(wp), PUBLIC ::   rn_tide_ramp_dt     !:
-   REAL(wp), PUBLIC ::   rn_scal_load        !:
+   REAL(wp), PUBLIC ::   rn_load        !:
    CHARACTER(lc), PUBLIC ::   cn_tide_load   !: 
    REAL(wp)         ::   rn_tide_gamma       ! Tidal tilt factor
@@ -124,5 +124,5 @@
       NAMELIST/nam_tide/ln_tide, nn_tide_var, ln_tide_dia, ln_tide_pot, rn_tide_gamma, &
          &              ln_scal_load, ln_read_load, cn_tide_load,         &
-         &              ln_tide_ramp, rn_scal_load, rn_tide_ramp_dt,      &
+         &              ln_tide_ramp, rn_load, rn_tide_ramp_dt,      &
          &              sn_tide_cnames
       !!----------------------------------------------------------------------
@@ -153,9 +153,9 @@
             WRITE(numout,*) '         Read load potential from file           ln_read_load    = ', ln_read_load
             WRITE(numout,*) '         Apply ramp on tides at startup          ln_tide_ramp    = ', ln_tide_ramp
-            WRITE(numout,*) '         Fraction of SSH used in scal. approx.   rn_scal_load    = ', rn_scal_load
+            WRITE(numout,*) '         Fraction of SSH used in scal. approx.   rn_load    = ', rn_load
             WRITE(numout,*) '         Duration (days) of ramp                 rn_tide_ramp_dt = ', rn_tide_ramp_dt
          ENDIF
       ELSE
-         rn_scal_load = 0._wp 
+         rn_load = 0._wp 
 
          IF(lwp) WRITE(numout,*)
@@ -193,5 +193,5 @@
       IF( nb_harmo == 0 )   CALL ctl_stop( 'tide_init : No tidal components set in nam_tide' )
       !
-      IF (.NOT.ln_scal_load ) rn_scal_load = 0._wp
+      IF (.NOT.ln_scal_load ) rn_load = 0._wp
       !
       ALLOCATE( amp_pot(jpi,jpj,nb_harmo),                      &
@@ -754,5 +754,5 @@
       !
       IF( ln_tide_dia ) THEN          ! Output total tidal potential (incl. load potential)
-         IF ( iom_use( "tide_pot" ) ) CALL iom_put( "tide_pot", pot_astro(:,:) + rn_scal_load * ssh(:,:,Kmm) )
+         IF ( iom_use( "tide_pot" ) ) CALL iom_put( "tide_pot", pot_astro(:,:) + rn_load * ssh(:,:,Kmm) )
       END IF
       !
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRA/eosbn2.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRA/eosbn2.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRA/eosbn2.F90	(revision 12443)
@@ -191,9 +191,9 @@
       !!                   ***  ROUTINE eos_insitu  ***
       !!
-      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from
+      !! ** Purpose :   Compute the in situ density (ratio rho/rho0) from
       !!       potential temperature and salinity using an equation of state
       !!       selected in the nameos namelist
       !!
-      !! ** Method  :   prd(t,s,z) = ( rho(t,s,z) - rau0 ) / rau0
+      !! ** Method  :   prd(t,s,z) = ( rho(t,s,z) - rho0 ) / rho0
       !!         with   prd    in situ density anomaly      no units
       !!                t      TEOS10: CT or EOS80: PT      Celsius
@@ -201,5 +201,5 @@
       !!                z      depth                        meters
       !!                rho    in situ density              kg/m^3
-      !!                rau0   reference density            kg/m^3
+      !!                rho0   reference density            kg/m^3
       !!
       !!     ln_teos10 : polynomial TEOS-10 equation of state is used for rho(t,s,z).
@@ -210,5 +210,5 @@
       !!
       !!     ln_seos : simplified equation of state
-      !!              prd(t,s,z) = ( -a0*(1+lambda/2*(T-T0)+mu*z+nu*(S-S0))*(T-T0) + b0*(S-S0) ) / rau0
+      !!              prd(t,s,z) = ( -a0*(1+lambda/2*(T-T0)+mu*z+nu*(S-S0))*(T-T0) + b0*(S-S0) ) / rho0
       !!              linear case function of T only: rn_alpha<>0, other coefficients = 0
       !!              linear eos function of T and S: rn_alpha and rn_beta<>0, other coefficients=0
@@ -267,5 +267,5 @@
             zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
             !
-            prd(ji,jj,jk) = (  zn * r1_rau0 - 1._wp  ) * ztm  ! density anomaly (masked)
+            prd(ji,jj,jk) = (  zn * r1_rho0 - 1._wp  ) * ztm  ! density anomaly (masked)
             !
          END_3D
@@ -283,5 +283,5 @@
                &  - rn_nu * zt * zs
                !                                 
-            prd(ji,jj,jk) = zn * r1_rau0 * ztm                ! density anomaly (masked)
+            prd(ji,jj,jk) = zn * r1_rho0 * ztm                ! density anomaly (masked)
          END_3D
          !
@@ -299,5 +299,5 @@
       !!                  ***  ROUTINE eos_insitu_pot  ***
       !!
-      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) and the
+      !! ** Purpose :   Compute the in situ density (ratio rho/rho0) and the
       !!      potential volumic mass (Kg/m3) from potential temperature and
       !!      salinity fields using an equation of state selected in the
@@ -379,5 +379,5 @@
                   prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp)                      ! potential density referenced at the surface
                   !
-                  prd(ji,jj,jk) = prd(ji,jj,jk) + (  zn_sto(jsmp) * r1_rau0 - 1._wp  )   ! density anomaly (masked)
+                  prd(ji,jj,jk) = prd(ji,jj,jk) + (  zn_sto(jsmp) * r1_rho0 - 1._wp  )   ! density anomaly (masked)
                END DO
                prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos
@@ -419,5 +419,5 @@
                prhop(ji,jj,jk) = zn0 * ztm                           ! potential density referenced at the surface
                !
-               prd(ji,jj,jk) = (  zn * r1_rau0 - 1._wp  ) * ztm      ! density anomaly (masked)
+               prd(ji,jj,jk) = (  zn * r1_rho0 - 1._wp  ) * ztm      ! density anomaly (masked)
             END_3D
          ENDIF
@@ -434,8 +434,8 @@
                &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs   &
                &  - rn_nu * zt * zs
-            prhop(ji,jj,jk) = ( rau0 + zn ) * ztm
+            prhop(ji,jj,jk) = ( rho0 + zn ) * ztm
             !                                                     ! density anomaly (masked)
             zn = zn - ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zh
-            prd(ji,jj,jk) = zn * r1_rau0 * ztm
+            prd(ji,jj,jk) = zn * r1_rho0 * ztm
             !
          END_3D
@@ -454,5 +454,5 @@
       !!                  ***  ROUTINE eos_insitu_2d  ***
       !!
-      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from
+      !! ** Purpose :   Compute the in situ density (ratio rho/rho0) from
       !!      potential temperature and salinity using an equation of state
       !!      selected in the nameos namelist. * 2D field case
@@ -508,5 +508,5 @@
             zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
             !
-            prd(ji,jj) = zn * r1_rau0 - 1._wp               ! unmasked in situ density anomaly
+            prd(ji,jj) = zn * r1_rho0 - 1._wp               ! unmasked in situ density anomaly
             !
          END_2D
@@ -524,5 +524,5 @@
                &  - rn_nu * zt * zs
                !
-            prd(ji,jj) = zn * r1_rau0               ! unmasked in situ density anomaly
+            prd(ji,jj) = zn * r1_rho0               ! unmasked in situ density anomaly
             !
          END_2D
@@ -588,5 +588,5 @@
             zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
             !
-            pab(ji,jj,jk,jp_tem) = zn * r1_rau0 * ztm
+            pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm
             !
             ! beta
@@ -609,5 +609,5 @@
             zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
             !
-            pab(ji,jj,jk,jp_sal) = zn / zs * r1_rau0 * ztm
+            pab(ji,jj,jk,jp_sal) = zn / zs * r1_rho0 * ztm
             !
          END_3D
@@ -622,8 +622,8 @@
             !
             zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
-            pab(ji,jj,jk,jp_tem) = zn * r1_rau0 * ztm   ! alpha
+            pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm   ! alpha
             !
             zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt
-            pab(ji,jj,jk,jp_sal) = zn * r1_rau0 * ztm   ! beta
+            pab(ji,jj,jk,jp_sal) = zn * r1_rho0 * ztm   ! beta
             !
          END_3D
@@ -694,5 +694,5 @@
             zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
             !
-            pab(ji,jj,jp_tem) = zn * r1_rau0
+            pab(ji,jj,jp_tem) = zn * r1_rho0
             !
             ! beta
@@ -715,5 +715,5 @@
             zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
             !
-            pab(ji,jj,jp_sal) = zn / zs * r1_rau0
+            pab(ji,jj,jp_sal) = zn / zs * r1_rho0
             !
             !
@@ -729,8 +729,8 @@
             !
             zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
-            pab(ji,jj,jp_tem) = zn * r1_rau0   ! alpha
+            pab(ji,jj,jp_tem) = zn * r1_rho0   ! alpha
             !
             zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt
-            pab(ji,jj,jp_sal) = zn * r1_rau0   ! beta
+            pab(ji,jj,jp_sal) = zn * r1_rho0   ! beta
             !
          END_2D
@@ -799,5 +799,5 @@
          zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
          !
-         pab(jp_tem) = zn * r1_rau0
+         pab(jp_tem) = zn * r1_rho0
          !
          ! beta
@@ -820,5 +820,5 @@
          zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
          !
-         pab(jp_sal) = zn / zs * r1_rau0
+         pab(jp_sal) = zn / zs * r1_rho0
          !
          !
@@ -831,8 +831,8 @@
          !
          zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
-         pab(jp_tem) = zn * r1_rau0   ! alpha
+         pab(jp_tem) = zn * r1_rho0   ! alpha
          !
          zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt
-         pab(jp_sal) = zn * r1_rau0   ! beta
+         pab(jp_sal) = zn * r1_rho0   ! beta
          !
       CASE DEFAULT
@@ -1052,10 +1052,10 @@
       !! ** Method  :   PE is defined analytically as the vertical 
       !!                   primitive of EOS times -g integrated between 0 and z>0.
-      !!                pen is the nonlinear bsq-PE anomaly: pen = ( PE - rau0 gz ) / rau0 gz - rd
+      !!                pen is the nonlinear bsq-PE anomaly: pen = ( PE - rho0 gz ) / rho0 gz - rd
       !!                                                      = 1/z * /int_0^z rd dz - rd 
       !!                                where rd is the density anomaly (see eos_rhd function)
       !!                ab_pe are partial derivatives of PE anomaly with respect to T and S:
-      !!                    ab_pe(1) = - 1/(rau0 gz) * dPE/dT + drd/dT = - d(pen)/dT
-      !!                    ab_pe(2) =   1/(rau0 gz) * dPE/dS + drd/dS =   d(pen)/dS
+      !!                    ab_pe(1) = - 1/(rho0 gz) * dPE/dT + drd/dT = - d(pen)/dT
+      !!                    ab_pe(2) =   1/(rho0 gz) * dPE/dS + drd/dS =   d(pen)/dS
       !!
       !! ** Action  : - pen         : PE anomaly given at T-points
@@ -1103,5 +1103,5 @@
             zn  = ( zn2 * zh + zn1 ) * zh + zn0
             !
-            ppen(ji,jj,jk)  = zn * zh * r1_rau0 * ztm
+            ppen(ji,jj,jk)  = zn * zh * r1_rho0 * ztm
             !
             ! alphaPE non-linear anomaly
@@ -1118,5 +1118,5 @@
             zn  = ( zn2 * zh + zn1 ) * zh + zn0
             !                              
-            pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rau0 * ztm
+            pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rho0 * ztm
             !
             ! betaPE non-linear anomaly
@@ -1133,5 +1133,5 @@
             zn  = ( zn2 * zh + zn1 ) * zh + zn0
             !                              
-            pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rau0 * ztm
+            pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rho0 * ztm
             !
          END_3D
@@ -1144,5 +1144,5 @@
             zh  = gdept(ji,jj,jk,Kmm)              ! depth in meters  at t-point
             ztm = tmask(ji,jj,jk)                ! tmask
-            zn  = 0.5_wp * zh * r1_rau0 * ztm
+            zn  = 0.5_wp * zh * r1_rho0 * ztm
             !                                    ! Potential Energy
             ppen(ji,jj,jk) = ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zn
@@ -1186,5 +1186,5 @@
       IF(lwm) WRITE( numond, nameos )
       !
-      rau0        = 1026._wp                 !: volumic mass of reference     [kg/m3]
+      rho0        = 1026._wp                 !: volumic mass of reference     [kg/m3]
       rcp         = 3991.86795711963_wp      !: heat capacity     [J/K]
       !
@@ -1598,5 +1598,5 @@
             WRITE(numout,*) '   ==>>>   use of simplified eos:    '
             WRITE(numout,*) '              rhd(dT=T-10,dS=S-35,Z) = [-a0*(1+lambda1/2*dT+mu1*Z)*dT '
-            WRITE(numout,*) '                                       + b0*(1+lambda2/2*dT+mu2*Z)*dS - nu*dT*dS] / rau0'
+            WRITE(numout,*) '                                       + b0*(1+lambda2/2*dT+mu2*Z)*dS - nu*dT*dS] / rho0'
             WRITE(numout,*) '              with the following coefficients :'
             WRITE(numout,*) '                 thermal exp. coef.    rn_a0      = ', rn_a0
@@ -1617,8 +1617,8 @@
       END SELECT
       !
-      rau0_rcp    = rau0 * rcp 
-      r1_rau0     = 1._wp / rau0
+      rho0_rcp    = rho0 * rcp 
+      r1_rho0     = 1._wp / rho0
       r1_rcp      = 1._wp / rcp
-      r1_rau0_rcp = 1._wp / rau0_rcp 
+      r1_rho0_rcp = 1._wp / rho0_rcp 
       !
       IF(lwp) THEN
@@ -1635,9 +1635,9 @@
       IF(lwp) WRITE(numout,*)
       IF(lwp) WRITE(numout,*) '   Associated physical constant'
-      IF(lwp) WRITE(numout,*) '      volumic mass of reference           rau0  = ', rau0   , ' kg/m^3'
-      IF(lwp) WRITE(numout,*) '      1. / rau0                        r1_rau0  = ', r1_rau0, ' m^3/kg'
+      IF(lwp) WRITE(numout,*) '      volumic mass of reference           rho0  = ', rho0   , ' kg/m^3'
+      IF(lwp) WRITE(numout,*) '      1. / rho0                        r1_rho0  = ', r1_rho0, ' m^3/kg'
       IF(lwp) WRITE(numout,*) '      ocean specific heat                 rcp   = ', rcp    , ' J/Kelvin'
-      IF(lwp) WRITE(numout,*) '      rau0 * rcp                       rau0_rcp = ', rau0_rcp
-      IF(lwp) WRITE(numout,*) '      1. / ( rau0 * rcp )           r1_rau0_rcp = ', r1_rau0_rcp
+      IF(lwp) WRITE(numout,*) '      rho0 * rcp                       rho0_rcp = ', rho0_rcp
+      IF(lwp) WRITE(numout,*) '      1. / ( rho0 * rcp )           r1_rho0_rcp = ', r1_rho0_rcp
       !
    END SUBROUTINE eos_init
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRA/traadv_fct.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRA/traadv_fct.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRA/traadv_fct.F90	(revision 12443)
@@ -20,5 +20,5 @@
    USE diaptr         ! poleward transport diagnostics
    USE diaar5         ! AR5 diagnostics
-   USE phycst  , ONLY : rau0_rcp
+   USE phycst  , ONLY : rho0_rcp
    USE zdf_oce , ONLY : ln_zad_Aimp
    !
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRA/traatf.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRA/traatf.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRA/traatf.F90	(revision 12443)
@@ -214,5 +214,5 @@
             ztd = pt(ji,jj,jk,jn,Kaa) - 2._wp * ztn + pt(ji,jj,jk,jn,Kbb)  ! time laplacian on tracers
             !
-            pt(ji,jj,jk,jn,Kmm) = ztn + atfp * ztd                      ! pt <-- filtered pt
+            pt(ji,jj,jk,jn,Kmm) = ztn + rn_atfp * ztd                      ! pt <-- filtered pt
          END_3D
          !
@@ -229,6 +229,6 @@
       !! 
       !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
-      !!             pt(Kmm)  = ( e3t(Kmm)*pt(Kmm) + atfp*[ e3t(Kbb)*pt(Kbb) - 2 e3t(Kmm)*pt(Kmm) + e3t_a*pt(Kaa) ] )
-      !!                       /( e3t(Kmm)         + atfp*[ e3t(Kbb)         - 2 e3t(Kmm)         + e3t(Kaa)    ] )
+      !!             pt(Kmm)  = ( e3t(Kmm)*pt(Kmm) + rn_atfp*[ e3t(Kbb)*pt(Kbb) - 2 e3t(Kmm)*pt(Kmm) + e3t_a*pt(Kaa) ] )
+      !!                       /( e3t(Kmm)         + rn_atfp*[ e3t(Kbb)         - 2 e3t(Kmm)         + e3t(Kaa)    ] )
       !!
       !! ** Action  : - pt(Kmm) ready for the next time step
@@ -272,6 +272,6 @@
       ENDIF
       zfact = 1._wp / p2dt
-      zfact1 = atfp * p2dt
-      zfact2 = zfact1 * r1_rau0
+      zfact1 = rn_atfp * p2dt
+      zfact2 = zfact1 * r1_rho0
       DO jn = 1, kjpt      
          DO_3D_00_00( 1, jpkm1 )
@@ -287,6 +287,6 @@
             ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
             !
-            ze3t_f = ze3t_n + atfp * ze3t_d
-            ztc_f  = ztc_n  + atfp * ztc_d
+            ze3t_f = ze3t_n + rn_atfp * ze3t_d
+            ztc_f  = ztc_n  + rn_atfp * ztc_d
             !
             ! Add asselin correction on scale factors:
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRA/trabbc.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRA/trabbc.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRA/trabbc.F90	(revision 12443)
@@ -66,5 +66,5 @@
       !!       ocean bottom can be computed once and is added to the temperature
       !!       trend juste above the bottom at each time step:
-      !!            ta = ta + Qsf / (rau0 rcp e3T) for k= mbkt
+      !!            ta = ta + Qsf / (rho0 rcp e3T) for k= mbkt
       !!       Where Qsf is the geothermal heat flux.
       !!
@@ -102,5 +102,5 @@
       ENDIF
       !
-      CALL iom_put ( "hfgeou" , rau0_rcp * qgh_trd0(:,:) )
+      CALL iom_put ( "hfgeou" , rho0_rcp * qgh_trd0(:,:) )
       IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' bbc  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
       !
@@ -162,5 +162,5 @@
          CASE ( 1 )                          !* constant flux
             IF(lwp) WRITE(numout,*) '   ==>>>   constant heat flux  =   ', rn_geoflx_cst
-            qgh_trd0(:,:) = r1_rau0_rcp * rn_geoflx_cst
+            qgh_trd0(:,:) = r1_rho0_rcp * rn_geoflx_cst
             !
          CASE ( 2 )                          !* variable geothermal heat flux : read the geothermal fluxes in mW/m2
@@ -179,5 +179,5 @@
 
             CALL fld_read( nit000, 1, sf_qgh )                         ! Read qgh data
-            qgh_trd0(:,:) = r1_rau0_rcp * sf_qgh(1)%fnow(:,:,1) * 1.e-3 ! conversion in W/m2
+            qgh_trd0(:,:) = r1_rho0_rcp * sf_qgh(1)%fnow(:,:,1) * 1.e-3 ! conversion in W/m2
             !
          CASE DEFAULT
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRA/tramle.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRA/tramle.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRA/tramle.F90	(revision 12443)
@@ -41,5 +41,5 @@
 
    REAL(wp) ::   r5_21 = 5.e0 / 21.e0   ! factor used in mle streamfunction computation
-   REAL(wp) ::   rb_c                   ! ML buoyancy criteria = g rho_c /rau0 where rho_c is defined in zdfmld
+   REAL(wp) ::   rb_c                   ! ML buoyancy criteria = g rho_c /rho0 where rho_c is defined in zdfmld
    REAL(wp) ::   rc_f                   ! MLE coefficient (= rn_ce / (5 km * fo) ) in nn_mle=1 case
 
@@ -112,5 +112,5 @@
          zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1  )  )    ! zc being 0 outside the ML t-points
          zmld(ji,jj) = zmld(ji,jj) + zc
-         zbm (ji,jj) = zbm (ji,jj) + zc * (rau0 - rhop(ji,jj,jk) ) * r1_rau0
+         zbm (ji,jj) = zbm (ji,jj) + zc * (rho0 - rhop(ji,jj,jk) ) * r1_rho0
          zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp
       END_3D
@@ -273,5 +273,5 @@
       IF( ln_mle ) THEN                ! MLE initialisation
          !
-         rb_c = grav * rn_rho_c_mle /rau0        ! Mixed Layer buoyancy criteria
+         rb_c = grav * rn_rho_c_mle /rho0        ! Mixed Layer buoyancy criteria
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) '      ML buoyancy criteria = ', rb_c, ' m/s2 '
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRA/traqsr.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRA/traqsr.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRA/traqsr.F90	(revision 12443)
@@ -87,5 +87,5 @@
       !!         I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) )
       !!         The temperature trend associated with the solar radiation penetration 
-      !!         is given by : zta = 1/e3t dk[ I ] / (rau0*Cp)
+      !!         is given by : zta = 1/e3t dk[ I ] / (rho0*Cp)
       !!         At the bottom, boudary condition for the radiation is no flux :
       !!      all heat which has not been absorbed in the above levels is put
@@ -155,5 +155,5 @@
          !
          DO jk = 1, nksr
-            qsr_hc(:,:,jk) = r1_rau0_rcp * ( etot3(:,:,jk) - etot3(:,:,jk+1) )
+            qsr_hc(:,:,jk) = r1_rho0_rcp * ( etot3(:,:,jk) - etot3(:,:,jk+1) )
          END DO
          !
@@ -227,5 +227,5 @@
          !
          DO_3D_00_00( 1, nksr )
-            qsr_hc(ji,jj,jk) = r1_rau0_rcp * ( zea(ji,jj,jk) - zea(ji,jj,jk+1) )
+            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zea(ji,jj,jk) - zea(ji,jj,jk+1) )
          END_3D
          !
@@ -234,6 +234,6 @@
       CASE( np_2BD  )            !==  2-bands fluxes  ==!
          !
-         zz0 =        rn_abs   * r1_rau0_rcp      ! surface equi-partition in 2-bands
-         zz1 = ( 1. - rn_abs ) * r1_rau0_rcp
+         zz0 =        rn_abs   * r1_rho0_rcp      ! surface equi-partition in 2-bands
+         zz1 = ( 1. - rn_abs ) * r1_rho0_rcp
          DO_3D_00_00( 1, nksr )
             zc0 = zz0 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi1r )
@@ -252,5 +252,5 @@
       ! sea-ice: store the 1st ocean level attenuation coefficient
       DO_2D_00_00
-         IF( qsr(ji,jj) /= 0._wp ) THEN   ;   fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) )
+         IF( qsr(ji,jj) /= 0._wp ) THEN   ;   fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rho0_rcp * qsr(ji,jj) )
          ELSE                             ;   fraqsr_1lev(ji,jj) = 1._wp
          ENDIF
@@ -262,5 +262,5 @@
          zetot(:,:,nksr+1:jpk) = 0._wp     ! below ~400m set to zero
          DO jk = nksr, 1, -1
-            zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) * rau0_rcp
+            zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) * rho0_rcp
          END DO         
          CALL iom_put( 'qsr3d', zetot )   ! 3D distribution of shortwave Radiation
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRA/trasbc.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRA/trasbc.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRA/trasbc.F90	(revision 12443)
@@ -124,11 +124,11 @@
       !                             !==  Now sbc tracer content fields  ==!
       DO_2D_01_00
-         sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj)   ! non solar heat flux
-         sbc_tsc(ji,jj,jp_sal) = r1_rau0     * sfx(ji,jj)   ! salt flux due to freezing/melting
+         sbc_tsc(ji,jj,jp_tem) = r1_rho0_rcp * qns(ji,jj)   ! non solar heat flux
+         sbc_tsc(ji,jj,jp_sal) = r1_rho0     * sfx(ji,jj)   ! salt flux due to freezing/melting
       END_2D
       IF( ln_linssh ) THEN                !* linear free surface  
          DO_2D_01_00
-            sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rau0 * emp(ji,jj) * pts(ji,jj,1,jp_tem,Kmm)
-            sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rau0 * emp(ji,jj) * pts(ji,jj,1,jp_sal,Kmm)
+            sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rho0 * emp(ji,jj) * pts(ji,jj,1,jp_tem,Kmm)
+            sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rho0 * emp(ji,jj) * pts(ji,jj,1,jp_sal,Kmm)
          END_2D
          IF( iom_use('emp_x_sst') )   CALL iom_put( "emp_x_sst", emp (:,:) * pts(:,:,1,jp_tem,Kmm) )
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRD/trddyn.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRD/trddyn.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRD/trddyn.F90	(revision 12443)
@@ -140,6 +140,6 @@
                               !                                    ! wind stress trends
                               ALLOCATE( z2dx(jpi,jpj) , z2dy(jpi,jpj) )
-                              z2dx(:,:) = ( utau_b(:,:) + utau(:,:) ) / ( e3u(:,:,1,Kmm) * rau0 )
-                              z2dy(:,:) = ( vtau_b(:,:) + vtau(:,:) ) / ( e3v(:,:,1,Kmm) * rau0 )
+                              z2dx(:,:) = ( utau_b(:,:) + utau(:,:) ) / ( e3u(:,:,1,Kmm) * rho0 )
+                              z2dy(:,:) = ( vtau_b(:,:) + vtau(:,:) ) / ( e3v(:,:,1,Kmm) * rho0 )
                               CALL iom_put( "utrd_tau", z2dx )
                               CALL iom_put( "vtrd_tau", z2dy )
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRD/trdglo.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRD/trdglo.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRD/trdglo.F90	(revision 12443)
@@ -76,5 +76,5 @@
       INTEGER ::   ji, jj, jk      ! dummy loop indices
       INTEGER ::   ikbu, ikbv      ! local integers
-      REAL(wp)::   zvm, zvt, zvs, z1_2rau0   ! local scalars
+      REAL(wp)::   zvm, zvt, zvs, z1_2rho0   ! local scalars
       REAL(wp), DIMENSION(jpi,jpj)  :: ztswu, ztswv, z2dx, z2dy   ! 2D workspace 
       !!----------------------------------------------------------------------
@@ -125,10 +125,10 @@
             !                 
             IF( ktrd == jpdyn_zdf ) THEN      ! zdf trend: compute separately the surface forcing trend
-               z1_2rau0 = 0.5_wp / rau0
+               z1_2rho0 = 0.5_wp / rho0
                DO_2D_10_10
                   zvt = ( utau_b(ji,jj) + utau(ji,jj) ) * tmask_i(ji+1,jj) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
-                     &                                                     * z1_2rau0       * e1e2u(ji,jj)
+                     &                                                     * z1_2rho0       * e1e2u(ji,jj)
                   zvs = ( vtau_b(ji,jj) + vtau(ji,jj) ) * tmask_i(ji,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
-                     &                                                     * z1_2rau0       * e1e2v(ji,jj)
+                     &                                                     * z1_2rho0       * e1e2v(ji,jj)
                   umo(jpdyn_tau) = umo(jpdyn_tau) + zvt
                   vmo(jpdyn_tau) = vmo(jpdyn_tau) + zvs
@@ -141,5 +141,5 @@
 !               !
 !               IF( ln_drgimp ) THEN                   ! implicit drag case: compute separately the bottom friction 
-!                  z1_2rau0 = 0.5_wp / rau0
+!                  z1_2rho0 = 0.5_wp / rho0
 !                  DO jj = 1, jpjm1
 !                     DO ji = 1, jpim1
@@ -203,5 +203,5 @@
          CALL eos( ts(:,:,:,:,Kmm), rhd, rhop )       ! now potential density
 
-         zcof = 0.5_wp / rau0             ! Density flux at w-point
+         zcof = 0.5_wp / rho0             ! Density flux at w-point
          zkz(:,:,1) = 0._wp
          DO jk = 2, jpk
@@ -209,5 +209,5 @@
          END DO
          
-         zcof   = 0.5_wp / rau0           ! Density flux at u and v-points
+         zcof   = 0.5_wp / rho0           ! Density flux at u and v-points
          DO_3D_10_10( 1, jpkm1 )
             zkx(ji,jj,jk) = zcof * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) * ( rhop(ji,jj,jk) + rhop(ji+1,jj,jk) )
@@ -347,5 +347,5 @@
  9546    FORMAT(' 0 < horizontal diffusion                                  : ', e20.13)
  9547    FORMAT(' 0 < vertical diffusion                                    : ', e20.13)
- 9548    FORMAT(' pressure gradient u2 = - 1/rau0 u.dz(rhop)                : ', e20.13, '  u.dz(rhop) =', e20.13)
+ 9548    FORMAT(' pressure gradient u2 = - 1/rho0 u.dz(rhop)                : ', e20.13, '  u.dz(rhop) =', e20.13)
          !
          ! Save potential to kinetic energy conversion for next time step
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRD/trdken.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRD/trdken.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRD/trdken.F90	(revision 12443)
@@ -102,5 +102,5 @@
       zke(:,1, : ) = 0._wp
       DO_3D_01_01( 1, jpkm1 )
-         zke(ji,jj,jk) = 0.5_wp * rau0 *( uu(ji  ,jj,jk,Kmm) * putrd(ji  ,jj,jk) * bu(ji  ,jj,jk)  &
+         zke(ji,jj,jk) = 0.5_wp * rho0 *( uu(ji  ,jj,jk,Kmm) * putrd(ji  ,jj,jk) * bu(ji  ,jj,jk)  &
             &                           + uu(ji-1,jj,jk,Kmm) * putrd(ji-1,jj,jk) * bu(ji-1,jj,jk)  &
             &                           + vv(ji,jj  ,jk,Kmm) * pvtrd(ji,jj  ,jk) * bv(ji,jj  ,jk)  &
@@ -123,5 +123,5 @@
                            zke2d(1,:) = 0._wp   ;   zke2d(:,1) = 0._wp
                            DO_2D_01_01
-                              zke2d(ji,jj) = r1_rau0 * 0.5_wp * (   z2dx(ji,jj) + z2dx(ji-1,jj)   &
+                              zke2d(ji,jj) = r1_rho0 * 0.5_wp * (   z2dx(ji,jj) + z2dx(ji-1,jj)   &
                               &                                   + z2dy(ji,jj) + z2dy(ji,jj-1)   ) * r1_bt(ji,jj,1)
                            END_2D
@@ -207,5 +207,5 @@
       !
       ! Local constant initialization 
-      zcoef = - rau0 * grav * 0.5_wp      
+      zcoef = - rho0 * grav * 0.5_wp      
       
       !  Surface value (also valid in partial step case)
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRD/trdvor.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRD/trdvor.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRD/trdvor.F90	(revision 12443)
@@ -105,6 +105,6 @@
          ztswu(:,:) = 0.e0   ;   ztswv(:,:) = 0.e0
          DO_2D_00_00
-            ztswu(ji,jj) = 0.5 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( e3u(ji,jj,1,Kmm) * rau0 )
-            ztswv(ji,jj) = 0.5 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( e3v(ji,jj,1,Kmm) * rau0 )
+            ztswu(ji,jj) = 0.5 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( e3u(ji,jj,1,Kmm) * rho0 )
+            ztswv(ji,jj) = 0.5 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( e3v(ji,jj,1,Kmm) * rho0 )
          END_2D
          !
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ZDF/zdfgls.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ZDF/zdfgls.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ZDF/zdfgls.F90	(revision 12443)
@@ -170,5 +170,5 @@
          !
          ! surface friction
-         ustar2_surf(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1)
+         ustar2_surf(ji,jj) = r1_rho0 * taum(ji,jj) * tmask(ji,jj,1)
          !   
 !!gm Rq we may add here r_ke0(_top/_bot) ?  ==>> think about that...
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ZDF/zdfiwm.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ZDF/zdfiwm.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ZDF/zdfiwm.F90	(revision 12443)
@@ -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 )
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ZDF/zdfmxl.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ZDF/zdfmxl.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ZDF/zdfmxl.F90	(revision 12443)
@@ -97,5 +97,5 @@
       nmln(:,:)  = nlb10               ! Initialization to the number of w ocean point
       hmlp(:,:)  = 0._wp               ! here hmlp used as a dummy variable, integrating vertically N^2
-      zN2_c = grav * rho_c * r1_rau0   ! convert density criteria into N^2 criteria
+      zN2_c = grav * rho_c * r1_rho0   ! convert density criteria into N^2 criteria
       DO_3D_11_11( nlb10, jpkm1 )
          ikt = mbkt(ji,jj)
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ZDF/zdfosm.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ZDF/zdfosm.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ZDF/zdfosm.F90	(revision 12443)
@@ -300,5 +300,5 @@
      DO_2D_00_00
         ! Surface downward irradiance (so always +ve)
-        zrad0(ji,jj) = qsr(ji,jj) * r1_rau0_rcp
+        zrad0(ji,jj) = qsr(ji,jj) * r1_rho0_rcp
         ! Downwards irradiance at base of boundary layer
         zradh(ji,jj) = zrad0(ji,jj) * ( zz0 * EXP( -hbl(ji,jj)/rn_si0 ) + zz1 * EXP( -hbl(ji,jj)/rn_si1) )
@@ -312,7 +312,7 @@
         zbeta    = rab_n(ji,jj,1,jp_sal)
         ! Upwards surface Temperature flux for non-local term
-        zwth0(ji,jj) = - qns(ji,jj) * r1_rau0_rcp * tmask(ji,jj,1)
+        zwth0(ji,jj) = - qns(ji,jj) * r1_rho0_rcp * tmask(ji,jj,1)
         ! Upwards surface salinity flux for non-local term
-        zws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * ts(ji,jj,1,jp_sal,Kmm)  + sfx(ji,jj) ) * r1_rau0 * tmask(ji,jj,1)
+        zws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * ts(ji,jj,1,jp_sal,Kmm)  + sfx(ji,jj) ) * r1_rho0 * tmask(ji,jj,1)
         ! Non radiative upwards surface buoyancy flux
         zwb0(ji,jj) = grav * zthermal * zwth0(ji,jj) -  grav * zbeta * zws0(ji,jj)
@@ -324,6 +324,6 @@
         zwbav(ji,jj) = grav  * zthermal * zwthav(ji,jj) - grav  * zbeta * zwsav(ji,jj)
         ! Surface upward velocity fluxes
-        zuw0(ji,jj) = -utau(ji,jj) * r1_rau0 * tmask(ji,jj,1)
-        zvw0(ji,jj) = -vtau(ji,jj) * r1_rau0 * tmask(ji,jj,1)
+        zuw0(ji,jj) = -utau(ji,jj) * r1_rho0 * tmask(ji,jj,1)
+        zvw0(ji,jj) = -vtau(ji,jj) * r1_rho0 * tmask(ji,jj,1)
         ! Friction velocity (zustar), at T-point : LMD94 eq. 2
         zustar(ji,jj) = MAX( SQRT( SQRT( zuw0(ji,jj) * zuw0(ji,jj) + zvw0(ji,jj) * zvw0(ji,jj) ) ), 1.0e-8 )
@@ -1250,10 +1250,10 @@
             IF ( iom_use("us_x") ) CALL iom_put( "us_x", tmask(:,:,1)*zustke*zcos_wind )   ! x surface Stokes drift
             IF ( iom_use("us_y") ) CALL iom_put( "us_y", tmask(:,:,1)*zustke*zsin_wind )  ! y surface Stokes drift
-            IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rau0*tmask(:,:,1)*zustar**2*zustke )
+            IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rho0*tmask(:,:,1)*zustar**2*zustke )
          ! Stokes drift read in from sbcwave  (=2).
          CASE(2)
             IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd )               ! x surface Stokes drift
             IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd )               ! y surface Stokes drift
-            IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rau0*tmask(:,:,1)*zustar**2* &
+            IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rho0*tmask(:,:,1)*zustar**2* &
                  & SQRT(ut0sd**2 + vt0sd**2 ) )
          END SELECT
@@ -1271,6 +1271,6 @@
          IF ( iom_use("zwstrl") ) CALL iom_put( "zwstrl", tmask(:,:,1)*zwstrl )         ! Langmuir velocity scale
          IF ( iom_use("zustar") ) CALL iom_put( "zustar", tmask(:,:,1)*zustar )         ! friction velocity scale
-         IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.*rau0*tmask(:,:,1)*zustar**3 ) ! BL depth internal to zdf_osm routine
-         IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*rau0*tmask(:,:,1)*zustar**2*zustke )
+         IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.*rho0*tmask(:,:,1)*zustar**3 ) ! BL depth internal to zdf_osm routine
+         IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*rho0*tmask(:,:,1)*zustar**2*zustke )
          IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl )               ! BL depth internal to zdf_osm routine
          IF ( iom_use("zhml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml )               ! ML depth internal to zdf_osm routine
@@ -1500,5 +1500,5 @@
      imld_rst(:,:)  = nlb10         ! Initialization to the number of w ocean point
      hbl(:,:)  = 0._wp              ! here hbl used as a dummy variable, integrating vertically N^2
-     zN2_c = grav * rho_c * r1_rau0 ! convert density criteria into N^2 criteria
+     zN2_c = grav * rho_c * r1_rho0 ! convert density criteria into N^2 criteria
      !
      hbl(:,:)  = 0._wp              ! here hbl used as a dummy variable, integrating vertically N^2
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ZDF/zdfric.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ZDF/zdfric.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ZDF/zdfric.F90	(revision 12443)
@@ -174,5 +174,5 @@
          !
          DO_2D_00_00
-            zustar = SQRT( taum(ji,jj) * r1_rau0 )
+            zustar = SQRT( taum(ji,jj) * r1_rho0 )
             zhek   = rn_ekmfc * zustar / ( ABS( ff_t(ji,jj) ) + rsmall )   ! Ekman depth
             zh_ekm(ji,jj) = MAX(  rn_mldmin , MIN( zhek , rn_mldmax )  )   ! set allowed range
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ZDF/zdftke.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ZDF/zdftke.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/ZDF/zdftke.F90	(revision 12443)
@@ -206,5 +206,5 @@
       !!--------------------------------------------------------------------
       !
-      zbbrau = rn_ebb / rau0       ! Local constant initialisation
+      zbbrau = rn_ebb / rho0       ! Local constant initialisation
       zfact1 = -.5_wp * rn_Dt 
       zfact2 = 1.5_wp * rn_Dt * rn_ediss
@@ -228,5 +228,5 @@
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !
-      !   en(bot)   = (ebb0/rau0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin)
+      !   en(bot)   = (ebb0/rho0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin)
       ! where ebb0 does not includes surface wave enhancement (i.e. ebb0=3.75)
       ! Note that stress averaged is done using an wet-only calculation of u and v at t-point like in zdfsh2
@@ -237,5 +237,5 @@
             zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
             zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )
-            !                       ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0)
+            !                       ! where 0.001875 = (rn_ebb0/rho0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0)
             zebot = - 0.001875_wp * rCdU_bot(ji,jj) * SQRT(  ( zmsku*( uu(ji,jj,mbkt(ji,jj),Kbb)+uu(ji-1,jj,mbkt(ji,jj),Kbb) ) )**2  &
                &                                           + ( zmskv*( vv(ji,jj,mbkt(ji,jj),Kbb)+vv(ji,jj-1,mbkt(ji,jj),Kbb) ) )**2  )
@@ -246,5 +246,5 @@
                zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) )
                zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) )
-               !                             ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000.  (CAUTION CdU<0)
+               !                             ! where 0.001875 = (rn_ebb0/rho0) * 0.5 = 3.75*0.5/1000.  (CAUTION CdU<0)
                zetop = - 0.001875_wp * rCdU_top(ji,jj) * SQRT(  ( zmsku*( uu(ji,jj,mikt(ji,jj),Kbb)+uu(ji-1,jj,mikt(ji,jj),Kbb) ) )**2  &
                   &                                           + ( zmskv*( vv(ji,jj,mikt(ji,jj),Kbb)+vv(ji,jj-1,mikt(ji,jj),Kbb) ) )**2  )
@@ -439,6 +439,6 @@
       zmxld(:,:,:)  = rmxl_min
       !
-      IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g)
-         zraug = vkarmn * 2.e5_wp / ( rau0 * grav )
+      IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rho0*g)
+         zraug = vkarmn * 2.e5_wp / ( rho0 * grav )
          DO_2D_00_00
             zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) )
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/oce.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/oce.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/oce.F90	(revision 12443)
@@ -28,5 +28,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)     ::   rn2b ,  rn2    !: brunt-vaisala frequency**2     [s-2]
    !
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   rhd    !: in situ density anomalie rhd=(rho-rau0)/rau0  [no units]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   rhd    !: in situ density anomalie rhd=(rho-rho0)/rho0  [no units]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   rhop   !: potential volumic mass                           [kg/m3]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   Cu_adv                   !: vertical Courant number (adaptive-implicit)
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OFF/dtadyn.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OFF/dtadyn.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OFF/dtadyn.F90	(revision 12443)
@@ -546,5 +546,5 @@
       ENDIF
 
-      ssh(:,:,Kmm) = ssh(:,:,Kmm) + atfp * ( ssh(:,:,Kbb) - 2 * ssh(:,:,Kmm) + ssh(:,:,Kaa))  
+      ssh(:,:,Kmm) = ssh(:,:,Kmm) + rn_atfp * ( ssh(:,:,Kbb) - 2 * ssh(:,:,Kmm) + ssh(:,:,Kaa))  
 
       !! Do we also need to time filter e3t??
@@ -629,5 +629,5 @@
       END DO
       !                                                ! Sea surface  elevation time-stepping
-      pssha(:,:) = ( psshb(:,:) - z2dt * ( r1_rau0 * pemp(:,:)  + zhdiv(:,:) ) ) * ssmask(:,:)
+      pssha(:,:) = ( psshb(:,:) - z2dt * ( r1_rho0 * pemp(:,:)  + zhdiv(:,:) ) ) * ssmask(:,:)
       !                                                 ! 
       !                                                 ! After acale factors at t-points ( z_star coordinate )
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/TOP/PISCES/P2Z/p2zexp.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/TOP/PISCES/P2Z/p2zexp.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/TOP/PISCES/P2Z/p2zexp.F90	(revision 12443)
@@ -122,5 +122,5 @@
         DO_2D_11_11
            zsedpocd = zsedpoca(ji,jj) - 2. * sedpocn(ji,jj) + sedpocb(ji,jj)      ! time laplacian on tracers
-           sedpocb(ji,jj) = sedpocn(ji,jj) + atfp * zsedpocd                     ! sedpocb <-- filtered sedpocn
+           sedpocb(ji,jj) = sedpocn(ji,jj) + rn_atfp * zsedpocd                     ! sedpocb <-- filtered sedpocn
            sedpocn(ji,jj) = zsedpoca(ji,jj)                                       ! sedpocn <-- sedpoca
         END_2D
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/TOP/TRP/trcatf.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/TOP/TRP/trcatf.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/TOP/TRP/trcatf.F90	(revision 12443)
@@ -71,5 +71,5 @@
       !!      the divergence of two consecutive time-steps and tr arrays
       !!      to prepare the next time_step:
-      !!         (tr(Kmm)) = (tr(Kmm)) + atfp [ (tr(Kbb)) + (tr(Kaa)) - 2 (tr(Kmm)) ]
+      !!         (tr(Kmm)) = (tr(Kmm)) + rn_atfp [ (tr(Kbb)) + (tr(Kaa)) - 2 (tr(Kmm)) ]
       !!
       !!
@@ -200,6 +200,6 @@
       !!                  /( e3t(:,:,:,Kmm)    + rbcp*[ e3t(:,:,:,Kbb)    - 2 e3t(:,:,:,Kmm)    + e3t(:,:,:,Kaa)    ] )   
       !!             ztm = 0                                                       otherwise
-      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
-      !!                  /( e3t(:,:,:,Kmm)    + atfp*[ e3t(:,:,:,Kbb)    - 2 e3t(:,:,:,Kmm)    + e3t(:,:,:,Kaa)    ] )
+      !!             tb  = ( e3t_n*tn + rn_atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
+      !!                  /( e3t(:,:,:,Kmm)    + rn_atfp*[ e3t(:,:,:,Kbb)    - 2 e3t(:,:,:,Kmm)    + e3t(:,:,:,Kaa)    ] )
       !!             tn  = ta 
       !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
@@ -222,6 +222,6 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
          IF( .NOT. ln_linssh ) THEN
-            rfact1 = atfp * rdttrc
-            rfact2 = rfact1 / rau0
+            rfact1 = rn_atfp * rdttrc
+            rfact2 = rfact1 / rho0
          ENDIF
         !  
@@ -241,6 +241,6 @@
             ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
             !
-            ze3t_f = ze3t_n + atfp * ze3t_d
-            ztc_f  = ztc_n  + atfp * ztc_d
+            ze3t_f = ze3t_n + rn_atfp * ze3t_d
+            ztc_f  = ztc_n  + rn_atfp * ztc_d
             !
             IF( .NOT. ln_linssh .AND. jk == mikt(ji,jj) ) THEN           ! first level 
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/TOP/TRP/trcsbc.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/TOP/TRP/trcsbc.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/TOP/TRP/trcsbc.F90	(revision 12443)
@@ -121,5 +121,5 @@
          DO jn = 1, jptra
             DO_2D_01_00
-               sbc_trc(ji,jj,jn) = zsfx(ji,jj) * r1_rau0 * ptr(ji,jj,1,jn,Kmm)
+               sbc_trc(ji,jj,jn) = zsfx(ji,jj) * r1_rho0 * ptr(ji,jj,1,jn,Kmm)
             END_2D
          END DO
@@ -129,5 +129,5 @@
          DO jn = 1, jptra
             DO_2D_01_00
-               sbc_trc(ji,jj,jn) = ( zsfx(ji,jj) + fmmflx(ji,jj) ) * r1_rau0 * ptr(ji,jj,1,jn,Kmm)
+               sbc_trc(ji,jj,jn) = ( zsfx(ji,jj) + fmmflx(ji,jj) ) * r1_rho0 * ptr(ji,jj,1,jn,Kmm)
             END_2D
          END DO
@@ -145,5 +145,5 @@
                ztfx  = zftra                        ! net tracer flux
                !
-               zdtra = r1_rau0 * ( ztfx + ( zsfx(ji,jj) + fmmflx(ji,jj) ) * ptr(ji,jj,1,jn,Kmm) ) 
+               zdtra = r1_rho0 * ( ztfx + ( zsfx(ji,jj) + fmmflx(ji,jj) ) * ptr(ji,jj,1,jn,Kmm) ) 
                IF ( zdtra < 0. ) THEN
                   zdtra  = MAX(zdtra, -ptr(ji,jj,1,jn,Kmm) * e3t(ji,jj,1,Kmm) / r2dttrc )   ! avoid negative concentrations to arise
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/TOP/oce_trc.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/TOP/oce_trc.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/TOP/oce_trc.F90	(revision 12443)
@@ -39,5 +39,5 @@
    USE oce , ONLY :   ts     =>    ts     !: 4D array contaning ( tn, sn )
    USE oce , ONLY :   rhop   =>    rhop   !: potential volumic mass (kg m-3) 
-   USE oce , ONLY :   rhd    =>    rhd    !: in situ density anomalie rhd=(rho-rau0)/rau0 (no units)
+   USE oce , ONLY :   rhd    =>    rhd    !: in situ density anomalie rhd=(rho-rho0)/rho0 (no units)
    USE oce , ONLY :   hdiv   =>    hdiv   !: horizontal divergence (1/s)
    USE oce , ONLY :   ssh    =>    ssh    !: sea surface height at t-point [m]   
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/TOP/trcbc.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/TOP/trcbc.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/TOP/trcbc.F90	(revision 12443)
@@ -416,5 +416,5 @@
             DO_2D_01_00
                DO jk = 1, nk_rnf(ji,jj)
-                  zrnf = (rnf(ji,jj) + rnf_b(ji,jj)) * 0.5_wp * r1_rau0 / h_rnf(ji,jj)
+                  zrnf = (rnf(ji,jj) + rnf_b(ji,jj)) * 0.5_wp * r1_rho0 / h_rnf(ji,jj)
                   ptr(ji,jj,jk,jn,Krhs) = ptr(ji,jj,jk,jn,Krhs)  + (ptr(ji,jj,jk,jn,Kmm) * zrnf)
                END DO
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/BENCH/EXPREF/namelist_cfg_orca025_like
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/BENCH/EXPREF/namelist_cfg_orca025_like	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/BENCH/EXPREF/namelist_cfg_orca025_like	(revision 12443)
@@ -189,5 +189,5 @@
    ln_dynspg_ts   = .true.   ! split-explicit free surface
       ln_bt_auto    = .false.    ! Number of sub-step defined from:
-         nn_baro      =  30         ! =F : the number of sub-step in rn_Dt seconds
+         nn_e      =  30         ! =F : the number of sub-step in rn_Dt seconds
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/BENCH/EXPREF/namelist_cfg_orca12_like
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/BENCH/EXPREF/namelist_cfg_orca12_like	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/BENCH/EXPREF/namelist_cfg_orca12_like	(revision 12443)
@@ -188,5 +188,5 @@
    ln_dynspg_ts   = .true.   ! split-explicit free surface
       ln_bt_auto    = .false.    ! Number of sub-step defined from:
-         nn_baro      =  30         ! =F : the number of sub-step in rn_Dt seconds
+         nn_e      =  30         ! =F : the number of sub-step in rn_Dt seconds
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/BENCH/EXPREF/namelist_cfg_orca1_like
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/BENCH/EXPREF/namelist_cfg_orca1_like	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/BENCH/EXPREF/namelist_cfg_orca1_like	(revision 12443)
@@ -188,5 +188,5 @@
    ln_dynspg_ts   = .true.   ! split-explicit free surface
       ln_bt_auto    = .false.    ! Number of sub-step defined from:
-         nn_baro      =  30         ! =F : the number of sub-step in rn_Dt seconds
+         nn_e      =  30         ! =F : the number of sub-step in rn_Dt seconds
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/BENCH/MY_SRC/diawri.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/BENCH/MY_SRC/diawri.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/BENCH/MY_SRC/diawri.F90	(revision 12443)
@@ -169,5 +169,5 @@
 
       IF ( iom_use("taubot") ) THEN                ! bottom stress
-         zztmp = rau0 * 0.25
+         zztmp = rho0 * 0.25
          z2d(:,:) = 0._wp
          DO jj = 2, jpjm1
@@ -214,5 +214,5 @@
       IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value
          ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
-         z2d(:,:) = rau0 * e1e2t(:,:)
+         z2d(:,:) = rho0 * e1e2t(:,:)
          DO jk = 1, jpk
             z3d(:,:,jk) = ww(:,:,jk) * z2d(:,:)
@@ -257,5 +257,5 @@
             END DO
          END DO
-         CALL iom_put( "heatc", rau0_rcp * z2d )   ! vertically integrated heat content (J/m2)
+         CALL iom_put( "heatc", rho0_rcp * z2d )   ! vertically integrated heat content (J/m2)
       ENDIF
 
@@ -269,5 +269,5 @@
             END DO
          END DO
-         CALL iom_put( "saltc", rau0 * z2d )          ! vertically integrated salt content (PSU*kg/m2)
+         CALL iom_put( "saltc", rho0 * z2d )          ! vertically integrated salt content (PSU*kg/m2)
       ENDIF
       !
@@ -295,5 +295,5 @@
          z2d(:,:) = 0.e0
          DO jk = 1, jpkm1
-            z3d(:,:,jk) = rau0 * uu(:,:,jk,Kmm) * e2u(:,:) * e3u(:,:,jk,Kmm) * umask(:,:,jk)
+            z3d(:,:,jk) = rho0 * uu(:,:,jk,Kmm) * e2u(:,:) * e3u(:,:,jk,Kmm) * umask(:,:,jk)
             z2d(:,:) = z2d(:,:) + z3d(:,:,jk)
          END DO
@@ -332,5 +332,5 @@
          z3d(:,:,jpk) = 0.e0
          DO jk = 1, jpkm1
-            z3d(:,:,jk) = rau0 * vv(:,:,jk,Kmm) * e1v(:,:) * e3v(:,:,jk,Kmm) * vmask(:,:,jk)
+            z3d(:,:,jk) = rho0 * vv(:,:,jk,Kmm) * e1v(:,:) * e3v(:,:,jk,Kmm) * vmask(:,:,jk)
          END DO
          CALL iom_put( "v_masstr", z3d )              ! mass transport in j-direction
@@ -373,5 +373,5 @@
          END DO
          CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
-         CALL iom_put( "tosmint", rau0 * z2d )        ! Vertical integral of temperature
+         CALL iom_put( "tosmint", rho0 * z2d )        ! Vertical integral of temperature
       ENDIF
       IF( iom_use("somint") ) THEN
@@ -385,5 +385,5 @@
          END DO
          CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
-         CALL iom_put( "somint", rau0 * z2d )         ! Vertical integral of salinity
+         CALL iom_put( "somint", rho0 * z2d )         ! Vertical integral of salinity
       ENDIF
 
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/BENCH/MY_SRC/zdfiwm.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/BENCH/MY_SRC/zdfiwm.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/BENCH/MY_SRC/zdfiwm.F90	(revision 12443)
@@ -85,10 +85,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
@@ -96,5 +96,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 
@@ -150,5 +150,5 @@
          DO ji = 1, jpi
             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 DO
@@ -179,5 +179,5 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               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 DO
          END DO
@@ -196,5 +196,5 @@
          DO jj= 1, jpj
             DO ji = 1, jpi
-               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 DO
          END DO
@@ -246,5 +246,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            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 DO
       END DO
@@ -259,5 +259,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)
@@ -305,5 +305,5 @@
          END DO
          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
@@ -349,5 +349,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) )
@@ -357,5 +357,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 )
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/CANAL/EXPREF/namelist_cfg
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/CANAL/EXPREF/namelist_cfg	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/CANAL/EXPREF/namelist_cfg	(revision 12443)
@@ -208,8 +208,8 @@
       ln_bt_av      = .true.     ! Time filtering of barotropic variables
          nn_bt_flt     = 1          ! Time filter choice  = 0 None
-         !                          !                     = 1 Boxcar over   nn_baro sub-steps
-         !                          !                     = 2 Boxcar over 2*nn_baro  "    "
+         !                          !                     = 1 Boxcar over   nn_e sub-steps
+         !                          !                     = 2 Boxcar over 2*nn_e  "    "
       ln_bt_auto    = .false.    ! Number of sub-step defined from:
-         nn_baro      =  24         ! =F : the number of sub-step in rn_Dt seconds
+         nn_e      =  24         ! =F : the number of sub-step in rn_Dt seconds
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/CANAL/MY_SRC/diawri.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/CANAL/MY_SRC/diawri.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/CANAL/MY_SRC/diawri.F90	(revision 12443)
@@ -169,5 +169,5 @@
 
       IF ( iom_use("taubot") ) THEN                ! bottom stress
-         zztmp = rau0 * 0.25
+         zztmp = rho0 * 0.25
          z2d(:,:) = 0._wp
          DO jj = 2, jpjm1
@@ -212,5 +212,5 @@
       IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value
          ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
-         z2d(:,:) = rau0 * e1e2t(:,:)
+         z2d(:,:) = rho0 * e1e2t(:,:)
          DO jk = 1, jpk
             z3d(:,:,jk) = ww(:,:,jk) * z2d(:,:)
@@ -272,5 +272,5 @@
             END DO
          END DO
-         CALL iom_put( "heatc", rau0_rcp * z2d )   ! vertically integrated heat content (J/m2)
+         CALL iom_put( "heatc", rho0_rcp * z2d )   ! vertically integrated heat content (J/m2)
       ENDIF
 
@@ -284,5 +284,5 @@
             END DO
          END DO
-         CALL iom_put( "saltc", rau0 * z2d )          ! vertically integrated salt content (PSU*kg/m2)
+         CALL iom_put( "saltc", rho0 * z2d )          ! vertically integrated salt content (PSU*kg/m2)
       ENDIF
       !
@@ -296,5 +296,5 @@
             END DO
          END DO
-         CALL iom_put( "salt2c", rau0 * z2d )          ! vertically integrated salt content (PSU*kg/m2)
+         CALL iom_put( "salt2c", rho0 * z2d )          ! vertically integrated salt content (PSU*kg/m2)
       ENDIF
       !
@@ -395,5 +395,5 @@
          z2d(:,:) = 0.e0
          DO jk = 1, jpkm1
-            z3d(:,:,jk) = rau0 * uu(:,:,jk,Kmm) * e2u(:,:) * e3u(:,:,jk,Kmm) * umask(:,:,jk)
+            z3d(:,:,jk) = rho0 * uu(:,:,jk,Kmm) * e2u(:,:) * e3u(:,:,jk,Kmm) * umask(:,:,jk)
             z2d(:,:) = z2d(:,:) + z3d(:,:,jk)
          END DO
@@ -432,5 +432,5 @@
          z3d(:,:,jpk) = 0.e0
          DO jk = 1, jpkm1
-            z3d(:,:,jk) = rau0 * vv(:,:,jk,Kmm) * e1v(:,:) * e3v(:,:,jk,Kmm) * vmask(:,:,jk)
+            z3d(:,:,jk) = rho0 * vv(:,:,jk,Kmm) * e1v(:,:) * e3v(:,:,jk,Kmm) * vmask(:,:,jk)
          END DO
          CALL iom_put( "v_masstr", z3d )              ! mass transport in j-direction
@@ -473,5 +473,5 @@
          END DO
          CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
-         CALL iom_put( "tosmint", rau0 * z2d )        ! Vertical integral of temperature
+         CALL iom_put( "tosmint", rho0 * z2d )        ! Vertical integral of temperature
       ENDIF
       IF( iom_use("somint") ) THEN
@@ -485,5 +485,5 @@
          END DO
          CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
-         CALL iom_put( "somint", rau0 * z2d )         ! Vertical integral of salinity
+         CALL iom_put( "somint", rho0 * z2d )         ! Vertical integral of salinity
       ENDIF
 
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/CANAL/MY_SRC/domvvl.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/CANAL/MY_SRC/domvvl.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/CANAL/MY_SRC/domvvl.F90	(revision 12443)
@@ -609,5 +609,5 @@
          ELSE
             tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 
-            &         + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
+            &         + rn_atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
          ENDIF
          tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/CANAL/MY_SRC/usrdef_istate.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/CANAL/MY_SRC/usrdef_istate.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/CANAL/MY_SRC/usrdef_istate.F90	(revision 12443)
@@ -218,5 +218,5 @@
          !
          zr_lambda2 = 1._wp / zlambda**2
-         zP0 = rau0 * zf0 * zumax * zlambda * SQRT(EXP(1._wp)/2._wp)
+         zP0 = rho0 * zf0 * zumax * zlambda * SQRT(EXP(1._wp)/2._wp)
          !
          DO jj=1, jpj
@@ -225,5 +225,5 @@
                zy = gphit(ji,jj) * 1.e3
                ! Surface pressure: P(x,y,z) = F(z) * Psurf(x,y)
-               zpsurf = zP0 * EXP(-(zx**2+zy**2)*zr_lambda2) - rau0 * ff_t(ji,jj) * rn_uzonal * zy
+               zpsurf = zP0 * EXP(-(zx**2+zy**2)*zr_lambda2) - rho0 * ff_t(ji,jj) * rn_uzonal * zy
                ! Sea level:
                pssh(ji,jj) = 0.
@@ -231,5 +231,5 @@
                   zdt = pssh(ji,jj)
                   zdzF = (1._wp - EXP(zdt-zH)) / (zH - 1._wp + EXP(-zH))   ! F'(z)
-                  zrho1 = rau0 * (1._wp + zn2*zdt/grav) - zdzF * zpsurf / grav    ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y)
+                  zrho1 = rho0 * (1._wp + zn2*zdt/grav) - zdzF * zpsurf / grav    ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y)
                   pssh(ji,jj) = zpsurf / (zrho1*grav) * ptmask(ji,jj,1)   ! ssh = Psurf / (Rho*g)
                END DO
@@ -237,11 +237,11 @@
                DO jk=1,jpk
                   zdt =  pdept(ji,jj,jk) 
-                  zrho1 = rau0 * (1._wp + zn2*zdt/grav)
+                  zrho1 = rho0 * (1._wp + zn2*zdt/grav)
                   IF (zdt < zH) THEN
                      zdzF = (1._wp-EXP(zdt-zH)) / (zH-1._wp + EXP(-zH))   ! F'(z)
                      zrho1 = zrho1 - zdzF * zpsurf / grav    ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y)
                   ENDIF
-                  !               pts(ji,jj,jk,jp_tem) = (20._wp + (rau0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk)
-                  pts(ji,jj,jk,jp_tem) = (10._wp + (rau0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk)
+                  !               pts(ji,jj,jk,jp_tem) = (20._wp + (rho0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk)
+                  pts(ji,jj,jk,jp_tem) = (10._wp + (rho0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk)
                END DO
             END DO
@@ -261,6 +261,6 @@
                   IF (zdu < zH) THEN
                      zf = (zH-1._wp-zdu+EXP(zdu-zH)) / (zH-1._wp+EXP(-zH))
-                     zdyPs = - za * zy * EXP(-(zx**2+zy**2)*zr_lambda2) - rau0 * ff_t(ji,jj) * rn_uzonal
-                     pu(ji,jj,jk) = - zf / ( rau0 * ff_t(ji,jj) ) * zdyPs * ptmask(ji,jj,jk) * ptmask(ji+1,jj,jk)
+                     zdyPs = - za * zy * EXP(-(zx**2+zy**2)*zr_lambda2) - rho0 * ff_t(ji,jj) * rn_uzonal
+                     pu(ji,jj,jk) = - zf / ( rho0 * ff_t(ji,jj) ) * zdyPs * ptmask(ji,jj,jk) * ptmask(ji+1,jj,jk)
                   ELSE
                      pu(ji,jj,jk) = 0._wp
@@ -279,5 +279,5 @@
                      zf = (zH-1._wp-zdv+EXP(zdv-zH)) / (zH-1._wp+EXP(-zH))
                      zdxPs = - za * zx * EXP(-(zx**2+zy**2)*zr_lambda2)
-                     pv(ji,jj,jk) = zf / ( rau0 * ff_f(ji,jj) ) * zdxPs * ptmask(ji,jj,jk) * ptmask(ji,jj+1,jk)
+                     pv(ji,jj,jk) = zf / ( rho0 * ff_f(ji,jj) ) * zdxPs * ptmask(ji,jj,jk) * ptmask(ji,jj+1,jk)
                   ELSE
                      pv(ji,jj,jk) = 0._wp
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/ISOMIP+/MY_SRC/eosbn2.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/ISOMIP+/MY_SRC/eosbn2.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/ISOMIP+/MY_SRC/eosbn2.F90	(revision 12443)
@@ -191,9 +191,9 @@
       !!                   ***  ROUTINE eos_insitu  ***
       !!
-      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from
+      !! ** Purpose :   Compute the in situ density (ratio rho/rho0) from
       !!       potential temperature and salinity using an equation of state
       !!       selected in the nameos namelist
       !!
-      !! ** Method  :   prd(t,s,z) = ( rho(t,s,z) - rau0 ) / rau0
+      !! ** Method  :   prd(t,s,z) = ( rho(t,s,z) - rho0 ) / rho0
       !!         with   prd    in situ density anomaly      no units
       !!                t      TEOS10: CT or EOS80: PT      Celsius
@@ -201,5 +201,5 @@
       !!                z      depth                        meters
       !!                rho    in situ density              kg/m^3
-      !!                rau0   reference density            kg/m^3
+      !!                rho0   reference density            kg/m^3
       !!
       !!     ln_teos10 : polynomial TEOS-10 equation of state is used for rho(t,s,z).
@@ -210,5 +210,5 @@
       !!
       !!     ln_seos : simplified equation of state
-      !!              prd(t,s,z) = ( -a0*(1+lambda/2*(T-T0)+mu*z+nu*(S-S0))*(T-T0) + b0*(S-S0) ) / rau0
+      !!              prd(t,s,z) = ( -a0*(1+lambda/2*(T-T0)+mu*z+nu*(S-S0))*(T-T0) + b0*(S-S0) ) / rho0
       !!              linear case function of T only: rn_alpha<>0, other coefficients = 0
       !!              linear eos function of T and S: rn_alpha and rn_beta<>0, other coefficients=0
@@ -216,5 +216,5 @@
       !!
       !!     ln_leos : linear ISOMIP equation of state
-      !!              prd(t,s,z) = ( -a0*(T-T0) + b0*(S-S0) ) / rau0
+      !!              prd(t,s,z) = ( -a0*(T-T0) + b0*(S-S0) ) / rho0
       !!              setup for ISOMIP linear eos
       !!
@@ -273,5 +273,5 @@
                   zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
                   !
-                  prd(ji,jj,jk) = (  zn * r1_rau0 - 1._wp  ) * ztm  ! density anomaly (masked)
+                  prd(ji,jj,jk) = (  zn * r1_rho0 - 1._wp  ) * ztm  ! density anomaly (masked)
                   !
                END DO
@@ -293,5 +293,5 @@
                      &  - rn_nu * zt * zs
                      !                                 
-                  prd(ji,jj,jk) = zn * r1_rau0 * ztm                ! density anomaly (masked)
+                  prd(ji,jj,jk) = zn * r1_rho0 * ztm                ! density anomaly (masked)
                END DO
             END DO
@@ -308,7 +308,7 @@
                   ztm = tmask(ji,jj,jk)
                   !
-                  zn =  rau0 * ( - rn_a0 * zt + rn_b0 * zs )
+                  zn =  rho0 * ( - rn_a0 * zt + rn_b0 * zs )
                   !                                 
-                  prd(ji,jj,jk) = zn * r1_rau0 * ztm                ! density anomaly (masked)
+                  prd(ji,jj,jk) = zn * r1_rho0 * ztm                ! density anomaly (masked)
                END DO
             END DO
@@ -328,5 +328,5 @@
       !!                  ***  ROUTINE eos_insitu_pot  ***
       !!
-      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) and the
+      !! ** Purpose :   Compute the in situ density (ratio rho/rho0) and the
       !!      potential volumic mass (Kg/m3) from potential temperature and
       !!      salinity fields using an equation of state selected in the
@@ -410,5 +410,5 @@
                         prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp)                      ! potential density referenced at the surface
                         !
-                        prd(ji,jj,jk) = prd(ji,jj,jk) + (  zn_sto(jsmp) * r1_rau0 - 1._wp  )   ! density anomaly (masked)
+                        prd(ji,jj,jk) = prd(ji,jj,jk) + (  zn_sto(jsmp) * r1_rho0 - 1._wp  )   ! density anomaly (masked)
                      END DO
                      prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos
@@ -454,5 +454,5 @@
                      prhop(ji,jj,jk) = zn0 * ztm                           ! potential density referenced at the surface
                      !
-                     prd(ji,jj,jk) = (  zn * r1_rau0 - 1._wp  ) * ztm      ! density anomaly (masked)
+                     prd(ji,jj,jk) = (  zn * r1_rho0 - 1._wp  ) * ztm      ! density anomaly (masked)
                   END DO
                END DO
@@ -473,8 +473,8 @@
                      &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs   &
                      &  - rn_nu * zt * zs
-                  prhop(ji,jj,jk) = ( rau0 + zn ) * ztm
+                  prhop(ji,jj,jk) = ( rho0 + zn ) * ztm
                   !                                                     ! density anomaly (masked)
                   zn = zn - ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zh
-                  prd(ji,jj,jk) = zn * r1_rau0 * ztm
+                  prd(ji,jj,jk) = zn * r1_rho0 * ztm
                   !
                END DO
@@ -492,8 +492,8 @@
                   ztm = tmask(ji,jj,jk)
                   !                                                     ! potential density referenced at the surface
-                  zn =  rau0 * ( - rn_a0 * zt + rn_b0 * zs )
-                  prhop(ji,jj,jk) = ( rau0 + zn ) * ztm
+                  zn =  rho0 * ( - rn_a0 * zt + rn_b0 * zs )
+                  prhop(ji,jj,jk) = ( rho0 + zn ) * ztm
                   !                                                     ! density anomaly (masked)
-                  prd(ji,jj,jk) = zn * r1_rau0 * ztm
+                  prd(ji,jj,jk) = zn * r1_rho0 * ztm
                   !
                END DO
@@ -514,5 +514,5 @@
       !!                  ***  ROUTINE eos_insitu_2d  ***
       !!
-      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from
+      !! ** Purpose :   Compute the in situ density (ratio rho/rho0) from
       !!      potential temperature and salinity using an equation of state
       !!      selected in the nameos namelist. * 2D field case
@@ -569,5 +569,5 @@
                zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
                !
-               prd(ji,jj) = zn * r1_rau0 - 1._wp               ! unmasked in situ density anomaly
+               prd(ji,jj) = zn * r1_rho0 - 1._wp               ! unmasked in situ density anomaly
                !
             END DO
@@ -589,5 +589,5 @@
                   &  - rn_nu * zt * zs
                   !
-               prd(ji,jj) = zn * r1_rau0               ! unmasked in situ density anomaly
+               prd(ji,jj) = zn * r1_rho0               ! unmasked in situ density anomaly
                !
             END DO
@@ -605,7 +605,7 @@
                zh    = pdep (ji,jj)                         ! depth at the partial step level
                !
-               zn =  rau0 * ( - rn_a0 * zt + rn_b0 * zs )
-                  !
-               prd(ji,jj) = zn * r1_rau0               ! unmasked in situ density anomaly
+               zn =  rho0 * ( - rn_a0 * zt + rn_b0 * zs )
+                  !
+               prd(ji,jj) = zn * r1_rho0               ! unmasked in situ density anomaly
                !
             END DO
@@ -676,5 +676,5 @@
                   zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
                   !
-                  pab(ji,jj,jk,jp_tem) = zn * r1_rau0 * ztm
+                  pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm
                   !
                   ! beta
@@ -697,5 +697,5 @@
                   zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
                   !
-                  pab(ji,jj,jk,jp_sal) = zn / zs * r1_rau0 * ztm
+                  pab(ji,jj,jk,jp_sal) = zn / zs * r1_rho0 * ztm
                   !
                END DO
@@ -714,8 +714,8 @@
                   !
                   zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
-                  pab(ji,jj,jk,jp_tem) = zn * r1_rau0 * ztm   ! alpha
+                  pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm   ! alpha
                   !
                   zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt
-                  pab(ji,jj,jk,jp_sal) = zn * r1_rau0 * ztm   ! beta
+                  pab(ji,jj,jk,jp_sal) = zn * r1_rho0 * ztm   ! beta
                   !
                END DO
@@ -733,9 +733,9 @@
                   ztm = tmask(ji,jj,jk)                   ! land/sea bottom mask = surf. mask
                   !
-                  zn  = rn_a0 * rau0
-                  pab(ji,jj,jk,jp_tem) = zn * r1_rau0 * ztm   ! alpha
-                  !
-                  zn  = rn_b0 * rau0
-                  pab(ji,jj,jk,jp_sal) = zn * r1_rau0 * ztm   ! beta
+                  zn  = rn_a0 * rho0
+                  pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm   ! alpha
+                  !
+                  zn  = rn_b0 * rho0
+                  pab(ji,jj,jk,jp_sal) = zn * r1_rho0 * ztm   ! beta
                   !
                END DO
@@ -809,5 +809,5 @@
                zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
                !
-               pab(ji,jj,jp_tem) = zn * r1_rau0
+               pab(ji,jj,jp_tem) = zn * r1_rho0
                !
                ! beta
@@ -830,5 +830,5 @@
                zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
                !
-               pab(ji,jj,jp_sal) = zn / zs * r1_rau0
+               pab(ji,jj,jp_sal) = zn / zs * r1_rho0
                !
                !
@@ -848,8 +848,8 @@
                !
                zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
-               pab(ji,jj,jp_tem) = zn * r1_rau0   ! alpha
+               pab(ji,jj,jp_tem) = zn * r1_rho0   ! alpha
                !
                zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt
-               pab(ji,jj,jp_sal) = zn * r1_rau0   ! beta
+               pab(ji,jj,jp_sal) = zn * r1_rho0   ! beta
                !
             END DO
@@ -867,9 +867,9 @@
                zh    = pdep (ji,jj)                   ! depth at the partial step level
                !
-               zn  = rn_a0 * rau0
-               pab(ji,jj,jp_tem) = zn * r1_rau0   ! alpha
-               !
-               zn  = rn_b0 * rau0
-               pab(ji,jj,jp_sal) = zn * r1_rau0   ! beta
+               zn  = rn_a0 * rho0
+               pab(ji,jj,jp_tem) = zn * r1_rho0   ! alpha
+               !
+               zn  = rn_b0 * rho0
+               pab(ji,jj,jp_sal) = zn * r1_rho0   ! beta
                !
             END DO
@@ -941,5 +941,5 @@
          zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
          !
-         pab(jp_tem) = zn * r1_rau0
+         pab(jp_tem) = zn * r1_rho0
          !
          ! beta
@@ -962,5 +962,5 @@
          zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
          !
-         pab(jp_sal) = zn / zs * r1_rau0
+         pab(jp_sal) = zn / zs * r1_rho0
          !
          !
@@ -973,8 +973,8 @@
          !
          zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
-         pab(jp_tem) = zn * r1_rau0   ! alpha
+         pab(jp_tem) = zn * r1_rho0   ! alpha
          !
          zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt
-         pab(jp_sal) = zn * r1_rau0   ! beta
+         pab(jp_sal) = zn * r1_rho0   ! beta
          !
       CASE( np_leos )                  !==  linear ISOMIP EOS  ==!
@@ -984,9 +984,9 @@
          zh    = pdep                    ! depth at the partial step level
          !
-         zn  = rn_a0 * rau0
-         pab(jp_tem) = zn * r1_rau0   ! alpha
-         !
-         zn  = rn_b0 * rau0
-         pab(jp_sal) = zn * r1_rau0   ! beta
+         zn  = rn_a0 * rho0
+         pab(jp_tem) = zn * r1_rho0   ! alpha
+         !
+         zn  = rn_b0 * rho0
+         pab(jp_sal) = zn * r1_rho0   ! beta
          !
       CASE DEFAULT
@@ -1214,10 +1214,10 @@
       !! ** Method  :   PE is defined analytically as the vertical 
       !!                   primitive of EOS times -g integrated between 0 and z>0.
-      !!                pen is the nonlinear bsq-PE anomaly: pen = ( PE - rau0 gz ) / rau0 gz - rd
+      !!                pen is the nonlinear bsq-PE anomaly: pen = ( PE - rho0 gz ) / rho0 gz - rd
       !!                                                      = 1/z * /int_0^z rd dz - rd 
       !!                                where rd is the density anomaly (see eos_rhd function)
       !!                ab_pe are partial derivatives of PE anomaly with respect to T and S:
-      !!                    ab_pe(1) = - 1/(rau0 gz) * dPE/dT + drd/dT = - d(pen)/dT
-      !!                    ab_pe(2) =   1/(rau0 gz) * dPE/dS + drd/dS =   d(pen)/dS
+      !!                    ab_pe(1) = - 1/(rho0 gz) * dPE/dT + drd/dT = - d(pen)/dT
+      !!                    ab_pe(2) =   1/(rho0 gz) * dPE/dS + drd/dS =   d(pen)/dS
       !!
       !! ** Action  : - pen         : PE anomaly given at T-points
@@ -1267,5 +1267,5 @@
                   zn  = ( zn2 * zh + zn1 ) * zh + zn0
                   !
-                  ppen(ji,jj,jk)  = zn * zh * r1_rau0 * ztm
+                  ppen(ji,jj,jk)  = zn * zh * r1_rho0 * ztm
                   !
                   ! alphaPE non-linear anomaly
@@ -1282,5 +1282,5 @@
                   zn  = ( zn2 * zh + zn1 ) * zh + zn0
                   !                              
-                  pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rau0 * ztm
+                  pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rho0 * ztm
                   !
                   ! betaPE non-linear anomaly
@@ -1297,5 +1297,5 @@
                   zn  = ( zn2 * zh + zn1 ) * zh + zn0
                   !                              
-                  pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rau0 * ztm
+                  pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rho0 * ztm
                   !
                END DO
@@ -1312,5 +1312,5 @@
                   zh  = gdept(ji,jj,jk,Kmm)              ! depth in meters  at t-point
                   ztm = tmask(ji,jj,jk)                ! tmask
-                  zn  = 0.5_wp * zh * r1_rau0 * ztm
+                  zn  = 0.5_wp * zh * r1_rho0 * ztm
                   !                                    ! Potential Energy
                   ppen(ji,jj,jk) = ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zn
@@ -1332,5 +1332,5 @@
                   zh  = gdept(ji,jj,jk,Kmm)                ! depth in meters  at t-point
                   ztm = tmask(ji,jj,jk)                  ! tmask
-                  zn  = 0.5_wp * zh * r1_rau0 * ztm
+                  zn  = 0.5_wp * zh * r1_rho0 * ztm
                   !                                    ! Potential Energy
                   ppen(ji,jj,jk) = 0.
@@ -1379,5 +1379,5 @@
       IF(lwm) WRITE( numond, nameos )
       !
-      rau0        = 1027.51_wp     !: volumic mass of reference     [kg/m3]
+      rho0        = 1027.51_wp     !: volumic mass of reference     [kg/m3]
       rcp         = 3974.00_wp     !: heat capacity     [J/K]
       !
@@ -1793,5 +1793,5 @@
             WRITE(numout,*) '   ==>>>   use of simplified eos:    '
             WRITE(numout,*) '              rhd(dT=T-10,dS=S-35,Z) = [-a0*(1+lambda1/2*dT+mu1*Z)*dT '
-            WRITE(numout,*) '                                       + b0*(1+lambda2/2*dT+mu2*Z)*dS - nu*dT*dS] / rau0'
+            WRITE(numout,*) '                                       + b0*(1+lambda2/2*dT+mu2*Z)*dS - nu*dT*dS] / rho0'
             WRITE(numout,*) '              with the following coefficients :'
             WRITE(numout,*) '                 thermal exp. coef.    rn_a0      = ', rn_a0
@@ -1810,5 +1810,5 @@
             WRITE(numout,*)
             WRITE(numout,*) '          use of linear ISOMIP eos:    rhd(dT=T-(-1),dS=S-(34.2),Z) = '
-            WRITE(numout,*) '             [ -a0*dT + b0*dS ]/rau0'
+            WRITE(numout,*) '             [ -a0*dT + b0*dS ]/rho0'
             WRITE(numout,*)
             WRITE(numout,*) '             thermal exp. coef.    rn_a0      = ', rn_a0
@@ -1822,8 +1822,8 @@
       END SELECT
       !
-      rau0_rcp    = rau0 * rcp 
-      r1_rau0     = 1._wp / rau0
+      rho0_rcp    = rho0 * rcp 
+      r1_rho0     = 1._wp / rho0
       r1_rcp      = 1._wp / rcp
-      r1_rau0_rcp = 1._wp / rau0_rcp 
+      r1_rho0_rcp = 1._wp / rho0_rcp 
       !
       IF(lwp) THEN
@@ -1840,9 +1840,9 @@
       IF(lwp) WRITE(numout,*)
       IF(lwp) WRITE(numout,*) '   Associated physical constant'
-      IF(lwp) WRITE(numout,*) '      volumic mass of reference           rau0  = ', rau0   , ' kg/m^3'
-      IF(lwp) WRITE(numout,*) '      1. / rau0                        r1_rau0  = ', r1_rau0, ' m^3/kg'
+      IF(lwp) WRITE(numout,*) '      volumic mass of reference           rho0  = ', rho0   , ' kg/m^3'
+      IF(lwp) WRITE(numout,*) '      1. / rho0                        r1_rho0  = ', r1_rho0, ' m^3/kg'
       IF(lwp) WRITE(numout,*) '      ocean specific heat                 rcp   = ', rcp    , ' J/Kelvin'
-      IF(lwp) WRITE(numout,*) '      rau0 * rcp                       rau0_rcp = ', rau0_rcp
-      IF(lwp) WRITE(numout,*) '      1. / ( rau0 * rcp )           r1_rau0_rcp = ', r1_rau0_rcp
+      IF(lwp) WRITE(numout,*) '      rho0 * rcp                       rho0_rcp = ', rho0_rcp
+      IF(lwp) WRITE(numout,*) '      1. / ( rho0 * rcp )           r1_rho0_rcp = ', r1_rho0_rcp
       !
    END SUBROUTINE eos_init
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/ISOMIP+/MY_SRC/sbcfwb.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/ISOMIP+/MY_SRC/sbcfwb.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/ISOMIP+/MY_SRC/sbcfwb.F90	(revision 12443)
@@ -122,5 +122,5 @@
             ! avoid the model to blow up for large ssh drop (isomip OCEAN3 with melt switch off and uniform T/S)
             IF (ln_isfcpl .AND. ln_isfcpl_cons) THEN
-               z_fwf = z_fwf + glob_sum( 'sbcfwb',  e1e2t(:,:) * risfcpl_cons_ssh(:,:) * rau0 )
+               z_fwf = z_fwf + glob_sum( 'sbcfwb',  e1e2t(:,:) * risfcpl_cons_ssh(:,:) * rho0 )
             END IF
             !
@@ -155,5 +155,5 @@
             a_fwb_b = a_fwb                           ! mean sea level taking into account the ice+snow
                                                       ! sum over the global domain
-            a_fwb   = glob_sum( 'sbcfwb', e1e2t(:,:) * ( ssh(:,:,Kmm) + snwice_mass(:,:) * r1_rau0 ) )
+            a_fwb   = glob_sum( 'sbcfwb', e1e2t(:,:) * ( ssh(:,:,Kmm) + snwice_mass(:,:) * r1_rho0 ) )
             a_fwb   = a_fwb * 1.e+3 / ( area * rday * 365. )     ! convert in Kg/m3/s = mm/s
 !!gm        !                                                      !!bug 365d year 
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT2_flux_cen2_cfg
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT2_flux_cen2_cfg	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT2_flux_cen2_cfg	(revision 12443)
@@ -144,8 +144,8 @@
       ln_bt_av      = .true.     ! Time filtering of barotropic variables
          nn_bt_flt     = 1          ! Time filter choice  = 0 None
-         !                          !                     = 1 Boxcar over   nn_baro sub-steps
-         !                          !                     = 2 Boxcar over 2*nn_baro  "    "
+         !                          !                     = 1 Boxcar over   nn_e sub-steps
+         !                          !                     = 2 Boxcar over 2*nn_e  "    "
       ln_bt_auto    = .true.     ! Number of sub-step defined from:
-         nn_baro      =  1         ! =F : the number of sub-step in rn_Dt seconds
+         nn_e      =  1         ! =F : the number of sub-step in rn_Dt seconds
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT2_flux_ubs_cfg
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT2_flux_ubs_cfg	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT2_flux_ubs_cfg	(revision 12443)
@@ -216,8 +216,8 @@
       ln_bt_av      = .true.     ! Time filtering of barotropic variables
          nn_bt_flt     = 1          ! Time filter choice  = 0 None
-         !                          !                     = 1 Boxcar over   nn_baro sub-steps
-         !                          !                     = 2 Boxcar over 2*nn_baro  "    "
+         !                          !                     = 1 Boxcar over   nn_e sub-steps
+         !                          !                     = 2 Boxcar over 2*nn_e  "    "
       ln_bt_auto    = .true.     ! Number of sub-step defined from:
-         nn_baro      =  1         ! =F : the number of sub-step in rn_Dt seconds
+         nn_e      =  1         ! =F : the number of sub-step in rn_Dt seconds
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT2_vect_eenH_cfg
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT2_vect_eenH_cfg	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT2_vect_eenH_cfg	(revision 12443)
@@ -144,8 +144,8 @@
       ln_bt_av      = .true.     ! Time filtering of barotropic variables
          nn_bt_flt     = 1          ! Time filter choice  = 0 None
-         !                          !                     = 1 Boxcar over   nn_baro sub-steps
-         !                          !                     = 2 Boxcar over 2*nn_baro  "    "
+         !                          !                     = 1 Boxcar over   nn_e sub-steps
+         !                          !                     = 2 Boxcar over 2*nn_e  "    "
       ln_bt_auto    = .true.     ! Number of sub-step defined from:
-         nn_baro      =  1         ! =F : the number of sub-step in rn_Dt seconds
+         nn_e      =  1         ! =F : the number of sub-step in rn_Dt seconds
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT2_vect_een_cfg
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT2_vect_een_cfg	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT2_vect_een_cfg	(revision 12443)
@@ -144,8 +144,8 @@
       ln_bt_av      = .true.     ! Time filtering of barotropic variables
          nn_bt_flt     = 1          ! Time filter choice  = 0 None
-         !                          !                     = 1 Boxcar over   nn_baro sub-steps
-         !                          !                     = 2 Boxcar over 2*nn_baro  "    "
+         !                          !                     = 1 Boxcar over   nn_e sub-steps
+         !                          !                     = 2 Boxcar over 2*nn_e  "    "
       ln_bt_auto    = .true.     ! Number of sub-step defined from:
-         nn_baro      =  1         ! =F : the number of sub-step in rn_Dt seconds
+         nn_e      =  1         ! =F : the number of sub-step in rn_Dt seconds
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT2_vect_ene_cfg
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT2_vect_ene_cfg	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT2_vect_ene_cfg	(revision 12443)
@@ -144,8 +144,8 @@
       ln_bt_av      = .true.     ! Time filtering of barotropic variables
          nn_bt_flt     = 1          ! Time filter choice  = 0 None
-         !                          !                     = 1 Boxcar over   nn_baro sub-steps
-         !                          !                     = 2 Boxcar over 2*nn_baro  "    "
+         !                          !                     = 1 Boxcar over   nn_e sub-steps
+         !                          !                     = 2 Boxcar over 2*nn_e  "    "
       ln_bt_auto    = .true.     ! Number of sub-step defined from:
-         nn_baro      =  1         ! =F : the number of sub-step in rn_Dt seconds
+         nn_e      =  1         ! =F : the number of sub-step in rn_Dt seconds
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT2_vect_ens_cfg
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT2_vect_ens_cfg	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT2_vect_ens_cfg	(revision 12443)
@@ -144,8 +144,8 @@
       ln_bt_av      = .true.     ! Time filtering of barotropic variables
          nn_bt_flt     = 1          ! Time filter choice  = 0 None
-         !                          !                     = 1 Boxcar over   nn_baro sub-steps
-         !                          !                     = 2 Boxcar over 2*nn_baro  "    "
+         !                          !                     = 1 Boxcar over   nn_e sub-steps
+         !                          !                     = 2 Boxcar over 2*nn_e  "    "
       ln_bt_auto    = .true.     ! Number of sub-step defined from:
-         nn_baro      =  1         ! =F : the number of sub-step in rn_Dt seconds
+         nn_e      =  1         ! =F : the number of sub-step in rn_Dt seconds
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT4_flux_cen2_cfg
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT4_flux_cen2_cfg	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT4_flux_cen2_cfg	(revision 12443)
@@ -144,8 +144,8 @@
       ln_bt_av      = .true.     ! Time filtering of barotropic variables
          nn_bt_flt     = 1          ! Time filter choice  = 0 None
-         !                          !                     = 1 Boxcar over   nn_baro sub-steps
-         !                          !                     = 2 Boxcar over 2*nn_baro  "    "
+         !                          !                     = 1 Boxcar over   nn_e sub-steps
+         !                          !                     = 2 Boxcar over 2*nn_e  "    "
       ln_bt_auto    = .true.     ! Number of sub-step defined from:
-         nn_baro      =  1         ! =F : the number of sub-step in rn_Dt seconds
+         nn_e      =  1         ! =F : the number of sub-step in rn_Dt seconds
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT4_flux_ubs_cfg
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT4_flux_ubs_cfg	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT4_flux_ubs_cfg	(revision 12443)
@@ -144,8 +144,8 @@
       ln_bt_av      = .true.     ! Time filtering of barotropic variables
          nn_bt_flt     = 1          ! Time filter choice  = 0 None
-         !                          !                     = 1 Boxcar over   nn_baro sub-steps
-         !                          !                     = 2 Boxcar over 2*nn_baro  "    "
+         !                          !                     = 1 Boxcar over   nn_e sub-steps
+         !                          !                     = 2 Boxcar over 2*nn_e  "    "
       ln_bt_auto    = .true.     ! Number of sub-step defined from:
-         nn_baro      =  1         ! =F : the number of sub-step in rn_Dt seconds
+         nn_e      =  1         ! =F : the number of sub-step in rn_Dt seconds
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT4_vect_eenH_cfg
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT4_vect_eenH_cfg	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT4_vect_eenH_cfg	(revision 12443)
@@ -144,8 +144,8 @@
       ln_bt_av      = .true.     ! Time filtering of barotropic variables
          nn_bt_flt     = 1          ! Time filter choice  = 0 None
-         !                          !                     = 1 Boxcar over   nn_baro sub-steps
-         !                          !                     = 2 Boxcar over 2*nn_baro  "    "
+         !                          !                     = 1 Boxcar over   nn_e sub-steps
+         !                          !                     = 2 Boxcar over 2*nn_e  "    "
       ln_bt_auto    = .true.     ! Number of sub-step defined from:
-         nn_baro      =  1         ! =F : the number of sub-step in rn_Dt seconds
+         nn_e      =  1         ! =F : the number of sub-step in rn_Dt seconds
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT4_vect_een_cfg
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT4_vect_een_cfg	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT4_vect_een_cfg	(revision 12443)
@@ -144,8 +144,8 @@
       ln_bt_av      = .true.     ! Time filtering of barotropic variables
          nn_bt_flt     = 1          ! Time filter choice  = 0 None
-         !                          !                     = 1 Boxcar over   nn_baro sub-steps
-         !                          !                     = 2 Boxcar over 2*nn_baro  "    "
+         !                          !                     = 1 Boxcar over   nn_e sub-steps
+         !                          !                     = 2 Boxcar over 2*nn_e  "    "
       ln_bt_auto    = .true.     ! Number of sub-step defined from:
-         nn_baro      =  1         ! =F : the number of sub-step in rn_Dt seconds
+         nn_e      =  1         ! =F : the number of sub-step in rn_Dt seconds
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT4_vect_ene_cfg
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT4_vect_ene_cfg	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT4_vect_ene_cfg	(revision 12443)
@@ -144,8 +144,8 @@
       ln_bt_av      = .true.     ! Time filtering of barotropic variables
          nn_bt_flt     = 1          ! Time filter choice  = 0 None
-         !                          !                     = 1 Boxcar over   nn_baro sub-steps
-         !                          !                     = 2 Boxcar over 2*nn_baro  "    "
+         !                          !                     = 1 Boxcar over   nn_e sub-steps
+         !                          !                     = 2 Boxcar over 2*nn_e  "    "
       ln_bt_auto    = .true.     ! Number of sub-step defined from:
-         nn_baro      =  1         ! =F : the number of sub-step in rn_Dt seconds
+         nn_e      =  1         ! =F : the number of sub-step in rn_Dt seconds
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT4_vect_ens_cfg
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT4_vect_ens_cfg	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/EXPREF/namelist_FCT4_vect_ens_cfg	(revision 12443)
@@ -144,8 +144,8 @@
       ln_bt_av      = .true.     ! Time filtering of barotropic variables
          nn_bt_flt     = 1          ! Time filter choice  = 0 None
-         !                          !                     = 1 Boxcar over   nn_baro sub-steps
-         !                          !                     = 2 Boxcar over 2*nn_baro  "    "
+         !                          !                     = 1 Boxcar over   nn_e sub-steps
+         !                          !                     = 2 Boxcar over 2*nn_e  "    "
       ln_bt_auto    = .true.     ! Number of sub-step defined from:
-         nn_baro      =  1         ! =F : the number of sub-step in rn_Dt seconds
+         nn_e      =  1         ! =F : the number of sub-step in rn_Dt seconds
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/MY_SRC/usrdef_istate.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/MY_SRC/usrdef_istate.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/LOCK_EXCHANGE/MY_SRC/usrdef_istate.F90	(revision 12443)
@@ -60,5 +60,5 @@
       !
       !  rn_a0 =  0.2   !  thermal expension coefficient (nn_eos= 1)
-      !  rho = rau0 - rn_a0 * (T-10) 
+      !  rho = rho0 - rn_a0 * (T-10) 
       !  delta_T = 25 degrees  ==>>  delta_rho = 25 * rn_a0 = 5 kg/m3
       !
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/OVERFLOW/EXPREF/namelist_sco_FCT2_flux_cen-ahm1000_cfg
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/OVERFLOW/EXPREF/namelist_sco_FCT2_flux_cen-ahm1000_cfg	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/OVERFLOW/EXPREF/namelist_sco_FCT2_flux_cen-ahm1000_cfg	(revision 12443)
@@ -155,8 +155,8 @@
       ln_bt_av      = .true.     ! Time filtering of barotropic variables
          nn_bt_flt     = 1          ! Time filter choice  = 0 None
-         !                          !                     = 1 Boxcar over   nn_baro sub-steps
-         !                          !                     = 2 Boxcar over 2*nn_baro  "    "
+         !                          !                     = 1 Boxcar over   nn_e sub-steps
+         !                          !                     = 2 Boxcar over 2*nn_e  "    "
       ln_bt_auto    = .true.     ! Number of sub-step defined from:
-         nn_baro      =  1         ! =F : the number of sub-step in rn_Dt seconds
+         nn_e      =  1         ! =F : the number of sub-step in rn_Dt seconds
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/OVERFLOW/EXPREF/namelist_sco_FCT2_flux_ubs_cfg
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/OVERFLOW/EXPREF/namelist_sco_FCT2_flux_ubs_cfg	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/OVERFLOW/EXPREF/namelist_sco_FCT2_flux_ubs_cfg	(revision 12443)
@@ -155,8 +155,8 @@
       ln_bt_av      = .true.     ! Time filtering of barotropic variables
          nn_bt_flt     = 1          ! Time filter choice  = 0 None
-         !                          !                     = 1 Boxcar over   nn_baro sub-steps
-         !                          !                     = 2 Boxcar over 2*nn_baro  "    "
+         !                          !                     = 1 Boxcar over   nn_e sub-steps
+         !                          !                     = 2 Boxcar over 2*nn_e  "    "
       ln_bt_auto    = .true.     ! Number of sub-step defined from:
-         nn_baro      =  1         ! =F : the number of sub-step in rn_Dt seconds
+         nn_e      =  1         ! =F : the number of sub-step in rn_Dt seconds
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/OVERFLOW/EXPREF/namelist_sco_FCT4_flux_cen-ahm1000_cfg
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/OVERFLOW/EXPREF/namelist_sco_FCT4_flux_cen-ahm1000_cfg	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/OVERFLOW/EXPREF/namelist_sco_FCT4_flux_cen-ahm1000_cfg	(revision 12443)
@@ -155,8 +155,8 @@
       ln_bt_av      = .true.     ! Time filtering of barotropic variables
          nn_bt_flt     = 1          ! Time filter choice  = 0 None
-         !                          !                     = 1 Boxcar over   nn_baro sub-steps
-         !                          !                     = 2 Boxcar over 2*nn_baro  "    "
+         !                          !                     = 1 Boxcar over   nn_e sub-steps
+         !                          !                     = 2 Boxcar over 2*nn_e  "    "
       ln_bt_auto    = .true.     ! Number of sub-step defined from:
-         nn_baro      =  1         ! =F : the number of sub-step in rn_Dt seconds
+         nn_e      =  1         ! =F : the number of sub-step in rn_Dt seconds
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/OVERFLOW/EXPREF/namelist_sco_FCT4_flux_ubs_cfg
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/OVERFLOW/EXPREF/namelist_sco_FCT4_flux_ubs_cfg	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/OVERFLOW/EXPREF/namelist_sco_FCT4_flux_ubs_cfg	(revision 12443)
@@ -155,8 +155,8 @@
       ln_bt_av      = .true.     ! Time filtering of barotropic variables
          nn_bt_flt     = 1          ! Time filter choice  = 0 None
-         !                          !                     = 1 Boxcar over   nn_baro sub-steps
-         !                          !                     = 2 Boxcar over 2*nn_baro  "    "
+         !                          !                     = 1 Boxcar over   nn_e sub-steps
+         !                          !                     = 2 Boxcar over 2*nn_e  "    "
       ln_bt_auto    = .true.     ! Number of sub-step defined from:
-         nn_baro      =  1         ! =F : the number of sub-step in rn_Dt seconds
+         nn_e      =  1         ! =F : the number of sub-step in rn_Dt seconds
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/OVERFLOW/EXPREF/namelist_zps_FCT2_flux_ubs_cfg
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/OVERFLOW/EXPREF/namelist_zps_FCT2_flux_ubs_cfg	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/OVERFLOW/EXPREF/namelist_zps_FCT2_flux_ubs_cfg	(revision 12443)
@@ -155,8 +155,8 @@
       ln_bt_av      = .true.     ! Time filtering of barotropic variables
          nn_bt_flt     = 1          ! Time filter choice  = 0 None
-         !                          !                     = 1 Boxcar over   nn_baro sub-steps
-         !                          !                     = 2 Boxcar over 2*nn_baro  "    "
+         !                          !                     = 1 Boxcar over   nn_e sub-steps
+         !                          !                     = 2 Boxcar over 2*nn_e  "    "
       ln_bt_auto    = .true.     ! Number of sub-step defined from:
-         nn_baro      =  1         ! =F : the number of sub-step in rn_Dt seconds
+         nn_e      =  1         ! =F : the number of sub-step in rn_Dt seconds
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/OVERFLOW/EXPREF/namelist_zps_FCT4_flux_ubs_cfg
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/OVERFLOW/EXPREF/namelist_zps_FCT4_flux_ubs_cfg	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/OVERFLOW/EXPREF/namelist_zps_FCT4_flux_ubs_cfg	(revision 12443)
@@ -216,8 +216,8 @@
       ln_bt_av      = .true.     ! Time filtering of barotropic variables
          nn_bt_flt     = 1          ! Time filter choice  = 0 None
-         !                          !                     = 1 Boxcar over   nn_baro sub-steps
-         !                          !                     = 2 Boxcar over 2*nn_baro  "    "
+         !                          !                     = 1 Boxcar over   nn_e sub-steps
+         !                          !                     = 2 Boxcar over 2*nn_e  "    "
       ln_bt_auto    = .true.     ! Number of sub-step defined from:
-         nn_baro      =  1         ! =F : the number of sub-step in rn_Dt seconds
+         nn_e      =  1         ! =F : the number of sub-step in rn_Dt seconds
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/OVERFLOW/EXPREF/namelist_zps_FCT4_vect_een_cfg
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/OVERFLOW/EXPREF/namelist_zps_FCT4_vect_een_cfg	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/OVERFLOW/EXPREF/namelist_zps_FCT4_vect_een_cfg	(revision 12443)
@@ -155,8 +155,8 @@
       ln_bt_av      = .true.     ! Time filtering of barotropic variables
          nn_bt_flt     = 1          ! Time filter choice  = 0 None
-         !                          !                     = 1 Boxcar over   nn_baro sub-steps
-         !                          !                     = 2 Boxcar over 2*nn_baro  "    "
+         !                          !                     = 1 Boxcar over   nn_e sub-steps
+         !                          !                     = 2 Boxcar over 2*nn_e  "    "
       ln_bt_auto    = .true.     ! Number of sub-step defined from:
-         nn_baro      =  1         ! =F : the number of sub-step in rn_Dt seconds
+         nn_e      =  1         ! =F : the number of sub-step in rn_Dt seconds
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/OVERFLOW/MY_SRC/usrdef_istate.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/OVERFLOW/MY_SRC/usrdef_istate.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/OVERFLOW/MY_SRC/usrdef_istate.F90	(revision 12443)
@@ -60,5 +60,5 @@
       !
       !  rn_a0 =  0.2   !  thermal expension coefficient (nn_eos= 1)
-      !  rho = rau0 - rn_a0 * (T-10) 
+      !  rho = rho0 - rn_a0 * (T-10) 
       !  delta_T = 10 degrees  ==>>  delta_rho = 10 * rn_a0 = 2 kg/m3
       !
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/VORTEX/EXPREF/1_namelist_cfg
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/VORTEX/EXPREF/1_namelist_cfg	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/VORTEX/EXPREF/1_namelist_cfg	(revision 12443)
@@ -213,8 +213,8 @@
       ln_bt_av      = .true.     ! Time filtering of barotropic variables
          nn_bt_flt     = 1          ! Time filter choice  = 0 None
-         !                          !                     = 1 Boxcar over   nn_baro sub-steps
-         !                          !                     = 2 Boxcar over 2*nn_baro  "    "
+         !                          !                     = 1 Boxcar over   nn_e sub-steps
+         !                          !                     = 2 Boxcar over 2*nn_e  "    "
       ln_bt_auto    = .false.    ! Number of sub-step defined from:
-         nn_baro      =  24         ! =F : the number of sub-step in rn_Dt seconds
+         nn_e      =  24         ! =F : the number of sub-step in rn_Dt seconds
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/VORTEX/EXPREF/namelist_cfg
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/VORTEX/EXPREF/namelist_cfg	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/VORTEX/EXPREF/namelist_cfg	(revision 12443)
@@ -204,8 +204,8 @@
       ln_bt_av      = .true.     ! Time filtering of barotropic variables
          nn_bt_flt     = 1          ! Time filter choice  = 0 None
-         !                          !                     = 1 Boxcar over   nn_baro sub-steps
-         !                          !                     = 2 Boxcar over 2*nn_baro  "    "
+         !                          !                     = 1 Boxcar over   nn_e sub-steps
+         !                          !                     = 2 Boxcar over 2*nn_e  "    "
       ln_bt_auto    = .false.    ! Number of sub-step defined from:
-         nn_baro      =  24         ! =F : the number of sub-step in rn_Dt seconds
+         nn_e      =  24         ! =F : the number of sub-step in rn_Dt seconds
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/VORTEX/MY_SRC/domvvl.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/VORTEX/MY_SRC/domvvl.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/VORTEX/MY_SRC/domvvl.F90	(revision 12443)
@@ -637,5 +637,5 @@
          ELSE
             tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 
-            &         + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
+            &         + rn_atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
          ENDIF
          tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/VORTEX/MY_SRC/usrdef_istate.F90
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/VORTEX/MY_SRC/usrdef_istate.F90	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/VORTEX/MY_SRC/usrdef_istate.F90	(revision 12443)
@@ -69,5 +69,5 @@
       zH = 0.5_wp * 5000._wp
       !
-      zP0 = rau0 * zf0 * zumax * zlambda * SQRT(EXP(1._wp)/2._wp)
+      zP0 = rho0 * zf0 * zumax * zlambda * SQRT(EXP(1._wp)/2._wp)
       !
       ! Sea level:
@@ -77,5 +77,5 @@
             zx = glamt(ji,jj) * 1.e3
             zy = gphit(ji,jj) * 1.e3
-            zrho1 = rau0 + za * EXP(-(zx**2+zy**2)/zlambda**2)
+            zrho1 = rho0 + za * EXP(-(zx**2+zy**2)/zlambda**2)
             pssh(ji,jj) = zP0 * EXP(-(zx**2+zy**2)/zlambda**2)/(zrho1*grav) * ptmask(ji,jj,1)
          END DO
@@ -89,10 +89,10 @@
             DO jk=1,jpk
                zdt =  pdept(ji,jj,jk) 
-               zrho1 = rau0 * (1._wp + zn2*zdt/grav)
+               zrho1 = rho0 * (1._wp + zn2*zdt/grav)
                IF (zdt < zH) THEN
                   zrho1 = zrho1 - zP0 * (1._wp-EXP(zdt-zH)) &
                           & * EXP(-(zx**2+zy**2)/zlambda**2) / (grav*(zH -1._wp + exp(-zH)));
                ENDIF
-               pts(ji,jj,jk,jp_tem) = (20._wp + (rau0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk)
+               pts(ji,jj,jk,jp_tem) = (20._wp + (rho0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk)
             END DO
          END DO
@@ -103,5 +103,5 @@
       !
       ! velocities:
-      za = 2._wp * zP0 / (zf0 * rau0 * zlambda**2)
+      za = 2._wp * zP0 / (zf0 * rho0 * zlambda**2)
       DO ji=1, jpim1
          DO jj=1, jpj
Index: /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/WAD/EXPREF/namelist_cfg
===================================================================
--- /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/WAD/EXPREF/namelist_cfg	(revision 12442)
+++ /NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/WAD/EXPREF/namelist_cfg	(revision 12443)
@@ -344,5 +344,5 @@
    ln_dynspg_ts  = .true.  !  split-explicit free surface
    ln_bt_auto    = .false. ! Number of sub-step defined from:
-   nn_baro       =  12     ! =F : the number of sub-step in rn_Dt seconds
+   nn_e       =  12     ! =F : the number of sub-step in rn_Dt seconds
 /
 !-----------------------------------------------------------------------
