Index: NEMO/trunk/src/OCE/TRA/eosbn2.F90
===================================================================
--- NEMO/trunk/src/OCE/TRA/eosbn2.F90	(revision 12377)
+++ NEMO/trunk/src/OCE/TRA/eosbn2.F90	(revision 12489)
@@ -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/trunk/src/OCE/TRA/traadv.F90
===================================================================
--- NEMO/trunk/src/OCE/TRA/traadv.F90	(revision 12377)
+++ NEMO/trunk/src/OCE/TRA/traadv.F90	(revision 12489)
@@ -92,9 +92,4 @@
       IF( ln_timing )   CALL timing_start('tra_adv')
       !
-      !                                          ! set time step
-      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt =         rdt   ! at nit000             (Euler)
-      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt = 2._wp * rdt   ! at nit000 or nit000+1 (Leapfrog)
-      ENDIF
-      !
       !                                         !==  effective transport  ==!
       zuu(:,:,jpk) = 0._wp
@@ -149,11 +144,11 @@
          CALL tra_adv_cen    ( kt, nit000, 'TRA',         zuu, zvv, zww, Kmm, pts, jpts, Krhs, nn_cen_h, nn_cen_v )
       CASE ( np_FCT )                                 ! FCT scheme      : 2nd / 4th order
-         CALL tra_adv_fct    ( kt, nit000, 'TRA', r2dt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_fct_h, nn_fct_v )
+         CALL tra_adv_fct    ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_fct_h, nn_fct_v )
       CASE ( np_MUS )                                 ! MUSCL
-         CALL tra_adv_mus    ( kt, nit000, 'TRA', r2dt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups ) 
+         CALL tra_adv_mus    ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups ) 
       CASE ( np_UBS )                                 ! UBS
-         CALL tra_adv_ubs    ( kt, nit000, 'TRA', r2dt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_ubs_v   )
+         CALL tra_adv_ubs    ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_ubs_v   )
       CASE ( np_QCK )                                 ! QUICKEST
-         CALL tra_adv_qck    ( kt, nit000, 'TRA', r2dt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs )
+         CALL tra_adv_qck    ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs )
       !
       END SELECT
Index: NEMO/trunk/src/OCE/TRA/traadv_fct.F90
===================================================================
--- NEMO/trunk/src/OCE/TRA/traadv_fct.F90	(revision 12377)
+++ NEMO/trunk/src/OCE/TRA/traadv_fct.F90	(revision 12489)
@@ -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/trunk/src/OCE/TRA/traatf.F90
===================================================================
--- NEMO/trunk/src/OCE/TRA/traatf.F90	(revision 12377)
+++ NEMO/trunk/src/OCE/TRA/traatf.F90	(revision 12489)
@@ -113,9 +113,4 @@
       IF( ln_bdy )   CALL bdy_tra( kt, Kbb, pts, Kaa )  ! BDY open boundaries
  
-      ! set time step size (Euler/Leapfrog)
-      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt =        rdt   ! at nit000             (Euler)
-      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt = 2._wp* rdt   ! at nit000 or nit000+1 (Leapfrog)
-      ENDIF
-
       ! trends computation initialisation
       IF( l_trdtra )   THEN                    
@@ -128,5 +123,5 @@
          ENDIF
          ! total trend for the non-time-filtered variables. 
-         zfact = 1.0 / rdt
+         zfact = 1.0 / rn_Dt
          ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from pts(Kmm) terms
          DO jk = 1, jpkm1
@@ -144,5 +139,5 @@
       ENDIF
 
-      IF( neuler == 0 .AND. kt == nit000 ) THEN       ! Euler time-stepping 
+      IF( l_1st_euler ) THEN       ! Euler time-stepping 
          !
          IF (l_trdtra .AND. .NOT. ln_linssh ) THEN   ! Zero Asselin filter contribution must be explicitly written out since for vvl
@@ -156,6 +151,6 @@
       ELSE                                            ! Leap-Frog + Asselin filter time stepping
          !
-         IF( ln_linssh ) THEN   ;   CALL tra_atf_fix( kt, Kbb, Kmm, Kaa, nit000,      'TRA', pts, jpts )  ! linear free surface 
-         ELSE                   ;   CALL tra_atf_vvl( kt, Kbb, Kmm, Kaa, nit000, rdt, 'TRA', pts, sbc_tsc, sbc_tsc_b, jpts )  ! non-linear free surface
+         IF( ln_linssh ) THEN   ;   CALL tra_atf_fix( kt, Kbb, Kmm, Kaa, nit000,        'TRA', pts, jpts )  ! linear free surface 
+         ELSE                   ;   CALL tra_atf_vvl( kt, Kbb, Kmm, Kaa, nit000, rn_Dt, 'TRA', pts, sbc_tsc, sbc_tsc_b, jpts )  ! non-linear free surface
          ENDIF
          !
@@ -167,5 +162,5 @@
       !
       IF( l_trdtra .AND. ln_linssh ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt     
-         zfact = 1._wp / r2dt             
+         zfact = 1._wp / rDt             
          DO jk = 1, jpkm1
             ztrdt(:,:,jk) = ( pts(:,:,jk,jp_tem,Kmm) - ztrdt(:,:,jk) ) * zfact
@@ -219,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
          !
@@ -234,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
@@ -277,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 )
@@ -292,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/trunk/src/OCE/TRA/trabbc.F90
===================================================================
--- NEMO/trunk/src/OCE/TRA/trabbc.F90	(revision 12377)
+++ NEMO/trunk/src/OCE/TRA/trabbc.F90	(revision 12489)
@@ -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/trunk/src/OCE/TRA/traldf_iso.F90
===================================================================
--- NEMO/trunk/src/OCE/TRA/traldf_iso.F90	(revision 12377)
+++ NEMO/trunk/src/OCE/TRA/traldf_iso.F90	(revision 12489)
@@ -109,5 +109,5 @@
       REAL(wp) ::  zmsku, zahu_w, zabe1, zcof1, zcoef3   ! local scalars
       REAL(wp) ::  zmskv, zahv_w, zabe2, zcof2, zcoef4   !   -      -
-      REAL(wp) ::  zcoef0, ze3w_2, zsign, z2dt, z1_2dt   !   -      -
+      REAL(wp) ::  zcoef0, ze3w_2, zsign                 !   -      -
       REAL(wp), DIMENSION(jpi,jpj)     ::   zdkt, zdk1t, z2d
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdit, zdjt, zftu, zftv, ztfw 
@@ -129,9 +129,4 @@
          &                        iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )   l_hst = .TRUE.
       !
-      !                                            ! set time step size (Euler/Leapfrog)
-      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   z2dt =     rdt      ! at nit000   (Euler)
-      ELSE                                        ;   z2dt = 2.* rdt      !             (Leapfrog)
-      ENDIF
-      z1_2dt = 1._wp / z2dt
       !
       IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0)
@@ -178,6 +173,6 @@
                DO_3D_10_10( 2, jpkm1 )
                   ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm)
-                  zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  )
-                  akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt
+                  zcoef0 = rDt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  )
+                  akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * r1_Dt
                END_3D
            ENDIF
Index: NEMO/trunk/src/OCE/TRA/traldf_triad.F90
===================================================================
--- NEMO/trunk/src/OCE/TRA/traldf_triad.F90	(revision 12377)
+++ NEMO/trunk/src/OCE/TRA/traldf_triad.F90	(revision 12489)
@@ -86,7 +86,7 @@
       INTEGER  ::  ip,jp,kp         ! dummy loop indices
       INTEGER  ::  ierr            ! local integer
-      REAL(wp) ::  zmsku, zabe1, zcof1, zcoef3          ! local scalars
-      REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4          !   -      -
-      REAL(wp) ::  zcoef0, ze3w_2, zsign, z2dt, z1_2dt  !   -      -
+      REAL(wp) ::  zmsku, zabe1, zcof1, zcoef3    ! local scalars
+      REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4    !   -      -
+      REAL(wp) ::  zcoef0, ze3w_2, zsign          !   -      -
       !
       REAL(wp) ::   zslope_skew, zslope_iso, zslope2, zbu, zbv
@@ -111,13 +111,9 @@
       l_hst = .FALSE.
       l_ptr = .FALSE.
-      IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) )      l_ptr = .TRUE. 
-      IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
-         &                        iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )   l_hst = .TRUE.
-      !
-      !                                                        ! set time step size (Euler/Leapfrog)
-      IF( neuler == 0 .AND. kt == kit000 ) THEN   ;   z2dt =     rdt      ! at nit000   (Euler)
-      ELSE                                        ;   z2dt = 2.* rdt      !             (Leapfrog)
+      IF( cdtype == 'TRA' ) THEN
+         IF( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf') )      l_ptr = .TRUE. 
+         IF( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR.                   &
+         &   iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  )   l_hst = .TRUE.
       ENDIF
-      z1_2dt = 1._wp / z2dt
       !
       IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0)
@@ -189,6 +185,6 @@
                DO_3D_10_10( 2, jpkm1 )
                   ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm)
-                  zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  )
-                  akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt
+                  zcoef0 = rDt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  )
+                  akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * r1_Dt
                END_3D
            ENDIF
Index: NEMO/trunk/src/OCE/TRA/tramle.F90
===================================================================
--- NEMO/trunk/src/OCE/TRA/tramle.F90	(revision 12377)
+++ NEMO/trunk/src/OCE/TRA/tramle.F90	(revision 12489)
@@ -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/trunk/src/OCE/TRA/tranpc.F90
===================================================================
--- NEMO/trunk/src/OCE/TRA/tranpc.F90	(revision 12377)
+++ NEMO/trunk/src/OCE/TRA/tranpc.F90	(revision 12489)
@@ -67,5 +67,5 @@
       LOGICAL  ::   l_bottom_reached, l_column_treated
       REAL(wp) ::   zta, zalfa, zsum_temp, zsum_alfa, zaw, zdz, zsum_z
-      REAL(wp) ::   zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw, z1_r2dt
+      REAL(wp) ::   zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw, z1_rDt
       REAL(wp), PARAMETER ::   zn2_zero = 1.e-14_wp             ! acceptance criteria for neutrality (N2==0)
       REAL(wp), DIMENSION(        jpk     )   ::   zvn2         ! vertical profile of N2 at 1 given point...
@@ -301,7 +301,7 @@
          !
          IF( l_trdtra ) THEN         ! send the Non penetrative mixing trends for diagnostic
-            z1_r2dt = 1._wp / (2._wp * rdt)
-            ztrdt(:,:,:) = ( pts(:,:,:,jp_tem,Kaa) - ztrdt(:,:,:) ) * z1_r2dt
-            ztrds(:,:,:) = ( pts(:,:,:,jp_sal,Kaa) - ztrds(:,:,:) ) * z1_r2dt
+            z1_rDt = 1._wp / (2._wp * rn_Dt)
+            ztrdt(:,:,:) = ( pts(:,:,:,jp_tem,Kaa) - ztrdt(:,:,:) ) * z1_rDt
+            ztrds(:,:,:) = ( pts(:,:,:,jp_sal,Kaa) - ztrds(:,:,:) ) * z1_rDt
             CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_npc, ztrdt )
             CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_npc, ztrds )
Index: NEMO/trunk/src/OCE/TRA/traqsr.F90
===================================================================
--- NEMO/trunk/src/OCE/TRA/traqsr.F90	(revision 12377)
+++ NEMO/trunk/src/OCE/TRA/traqsr.F90	(revision 12489)
@@ -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
@@ -135,6 +135,5 @@
       !                         !-----------------------------------!
       IF( kt == nit000 ) THEN          !==  1st time step  ==!
-!!gm case neuler  not taken into account....
-         IF( ln_rstart .AND. iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 ) THEN    ! read in restart
+         IF( ln_rstart .AND. iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0  .AND. .NOT.l_1st_euler ) THEN    ! read in restart
             IF(lwp) WRITE(numout,*) '          nit000-1 qsr tracer content forcing field read in the restart file'
             z1_2 = 0.5_wp
@@ -156,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
          !
@@ -228,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
          !
@@ -235,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 )
@@ -253,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
@@ -263,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/trunk/src/OCE/TRA/trasbc.F90
===================================================================
--- NEMO/trunk/src/OCE/TRA/trasbc.F90	(revision 12377)
+++ NEMO/trunk/src/OCE/TRA/trasbc.F90	(revision 12489)
@@ -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/trunk/src/OCE/TRA/trazdf.F90
===================================================================
--- NEMO/trunk/src/OCE/TRA/trazdf.F90	(revision 12377)
+++ NEMO/trunk/src/OCE/TRA/trazdf.F90	(revision 12489)
@@ -66,8 +66,4 @@
       ENDIF
       !
-      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt =      rdt   ! at nit000, =   rdt (restarting with Euler time stepping)
-      ELSEIF( kt <= nit000 + 1           ) THEN   ;   r2dt = 2. * rdt   ! otherwise, = 2 rdt (leapfrog)
-      ENDIF
-      !
       IF( l_trdtra )   THEN                  !* Save ta and sa trends
          ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
@@ -77,5 +73,5 @@
       !
       !                                      !* compute lateral mixing trend and add it to the general trend
-      CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt, Kbb, Kmm, Krhs, pts, Kaa, jpts ) 
+      CALL tra_zdf_imp( kt, nit000, 'TRA', rDt, Kbb, Kmm, Krhs, pts, Kaa, jpts ) 
 
 !!gm WHY here !   and I don't like that !
@@ -89,7 +85,7 @@
          DO jk = 1, jpkm1
             ztrdt(:,:,jk) = ( ( pts(:,:,jk,jp_tem,Kaa)*e3t(:,:,jk,Kaa) - pts(:,:,jk,jp_tem,Kbb)*e3t(:,:,jk,Kbb) ) &
-               &          / (e3t(:,:,jk,Kmm)*r2dt) ) - ztrdt(:,:,jk)
+               &          / (e3t(:,:,jk,Kmm)*rDt) ) - ztrdt(:,:,jk)
             ztrds(:,:,jk) = ( ( pts(:,:,jk,jp_sal,Kaa)*e3t(:,:,jk,Kaa) - pts(:,:,jk,jp_sal,Kbb)*e3t(:,:,jk,Kbb) ) &
-              &           / (e3t(:,:,jk,Kmm)*r2dt) ) - ztrds(:,:,jk)
+              &           / (e3t(:,:,jk,Kmm)*rDt) ) - ztrds(:,:,jk)
          END DO
 !!gm this should be moved in trdtra.F90 and done on all trends
