Index: /branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/tranxt.F90
===================================================================
--- /branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 1846)
+++ /branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 1847)
@@ -15,4 +15,5 @@
    !!            3.0  !  2008-06  (G. Madec)  time stepping always done in trazdf
    !!            3.1  !  2009-02  (G. Madec, R. Benshila)  re-introduce the vvl option
+   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  semi-implicit hpg with asselin filter
    !!----------------------------------------------------------------------
 
@@ -44,10 +45,11 @@
    PUBLIC   tra_nxt    ! routine called by step.F90
 
-   REAL(wp), DIMENSION(jpk) ::   r2dt_t   ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler)
+   REAL(wp)                 ::   rbcp, r1_2bcp   ! Brown & Campana parameters for semi-implicit hpg
+   REAL(wp), DIMENSION(jpk) ::   r2dt_t          ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler)
 
    !! * Substitutions
 #  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 
+   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -85,6 +87,6 @@
       INTEGER, INTENT(in) ::   kt    ! ocean time-step index
       !!
-      INTEGER  ::   jk    ! dummy loop indices
-      REAL(wp) ::   zfact ! temporary scalars
+      INTEGER  ::   jk      ! dummy loop indices
+      REAL(wp) ::   zfact   ! temporary scalars
       !!----------------------------------------------------------------------
 
@@ -93,4 +95,7 @@
          IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap'
          IF(lwp) WRITE(numout,*) '~~~~~~~'
+         !
+         rbcp    = 0.25 * (1 + atfp) * (1 + atfp * atfp)       ! Brown & Campana parameter for semi-implicit hpg
+         r1_2bcp = 1.e0 - 2.e0 * rbcp
       ENDIF
 
@@ -160,7 +165,8 @@
       !!             hydrostatic pressure gradient (ln_dynhpg_imp = T)
       !!              - swap tracer fields to prepare the next time_step.
-      !!                This can be summurized for tempearture as:
-      !!             ztm = (ta+2tn+tb)/4       ln_dynhpg_imp = T
-      !!             ztm = 0                   otherwise
+      !!                This can be summurized for temperature as:
+      !!             ztm = (rbcp*ta+(2-rbcp)*tn+rbcp*tb)       ln_dynhpg_imp = T
+      !!             ztm = 0                                   otherwise
+      !!                   with rbcp=(1+atfp)*(1+atfp*atfp)/4
       !!             tb  = tn + atfp*[ tb - 2 tn + ta ]
       !!             tn  = ta 
@@ -189,24 +195,20 @@
          IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step
             DO jk = 1, jpkm1                              ! (only swap)
+               tn(:,:,jk) = ta(:,:,jk)                                                          ! tb <-- tn
+               sn(:,:,jk) = sa(:,:,jk)
+            END DO
+         ELSE                                             ! general case (Leapfrog + Asselin filter)
+            DO jk = 1, jpkm1
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     tn(ji,jj,jk) = ta(ji,jj,jk)                                           ! tb <-- tn
+                     ztm = rbcp * ( ta(ji,jj,jk) + tb(ji,jj,jk) ) + r1_2bcp * tn(ji,jj,jk)      ! mean t
+                     zsm = rbcp * ( sa(ji,jj,jk) + sb(ji,jj,jk) ) + r1_2bcp * sn(ji,jj,jk) 
+                     ztf = atfp * ( ta(ji,jj,jk) + tb(ji,jj,jk)   - 2.      * tn(ji,jj,jk) )    ! Asselin filter on t 
+                     zsf = atfp * ( sa(ji,jj,jk) + sb(ji,jj,jk)   - 2.      * sn(ji,jj,jk) )
+                     tb(ji,jj,jk) = tn(ji,jj,jk) + ztf                                          ! tb <-- filtered tn 
+                     sb(ji,jj,jk) = sn(ji,jj,jk) + zsf
+                     tn(ji,jj,jk) = ta(ji,jj,jk)                                                ! tn <-- ta
                      sn(ji,jj,jk) = sa(ji,jj,jk)
-                  END DO
-               END DO
-            END DO
-         ELSE                                             ! general case (Leapfrog + Asselin filter
-            DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     ztm = 0.25 * ( ta(ji,jj,jk) + 2.* tn(ji,jj,jk) + tb(ji,jj,jk) )       ! mean t
-                     zsm = 0.25 * ( sa(ji,jj,jk) + 2.* sn(ji,jj,jk) + sb(ji,jj,jk) )
-                     ztf = atfp * ( ta(ji,jj,jk) - 2.* tn(ji,jj,jk) + tb(ji,jj,jk) )       ! Asselin filter on t 
-                     zsf = atfp * ( sa(ji,jj,jk) - 2.* sn(ji,jj,jk) + sb(ji,jj,jk) )
-                     tb(ji,jj,jk) = tn(ji,jj,jk) + ztf                                     ! tb <-- filtered tn 
-                     sb(ji,jj,jk) = sn(ji,jj,jk) + zsf
-                     tn(ji,jj,jk) = ta(ji,jj,jk)                                           ! tn <-- ta
-                     sn(ji,jj,jk) = sa(ji,jj,jk)
-                     ta(ji,jj,jk) = ztm                                                    ! ta <-- mean t
+                     ta(ji,jj,jk) = ztm                                                         ! ta <-- mean t
                      sa(ji,jj,jk) = zsm
                   END DO
@@ -220,20 +222,16 @@
          IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step
             DO jk = 1, jpkm1
+               tn(:,:,jk) = ta(:,:,jk)                                                    ! tn <-- ta
+               sn(:,:,jk) = sa(:,:,jk)
+            END DO
+         ELSE                                             ! general case (Leapfrog + Asselin filter
+            DO jk = 1, jpkm1
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     tn(ji,jj,jk) = ta(ji,jj,jk)                                           ! tn <-- ta
-                     sn(ji,jj,jk) = sa(ji,jj,jk)
-                  END DO
-               END DO
-            END DO
-         ELSE                                             ! general case (Leapfrog + Asselin filter
-            DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     ztf = atfp * ( ta(ji,jj,jk) - 2.* tn(ji,jj,jk) + tb(ji,jj,jk) )       ! Asselin filter on t 
+                     ztf = atfp * ( ta(ji,jj,jk) - 2.* tn(ji,jj,jk) + tb(ji,jj,jk) )      ! Asselin filter on t 
                      zsf = atfp * ( sa(ji,jj,jk) - 2.* sn(ji,jj,jk) + sb(ji,jj,jk) )
-                     tb(ji,jj,jk) = tn(ji,jj,jk) + ztf                                     ! tb <-- filtered tn 
+                     tb(ji,jj,jk) = tn(ji,jj,jk) + ztf                                    ! tb <-- filtered tn 
                      sb(ji,jj,jk) = sn(ji,jj,jk) + zsf
-                     tn(ji,jj,jk) = ta(ji,jj,jk)                                           ! tn <-- ta
+                     tn(ji,jj,jk) = ta(ji,jj,jk)                                          ! tn <-- ta
                      sn(ji,jj,jk) = sa(ji,jj,jk)
                   END DO
@@ -258,8 +256,8 @@
       !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
       !!              - swap tracer fields to prepare the next time_step.
-      !!                This can be summurized for tempearture as:
-      !!             ztm = (e3t_a*ta+2*e3t_n*tn+e3t_b*tb)   ln_dynhpg_imp = T
-      !!                  /(e3t_a   +2*e3t_n   +e3t_b   )    
-      !!             ztm = 0                                otherwise
+      !!                This can be summurized for temperature as:
+      !!             ztm = ( rbcp*e3t_a*ta + (2-rbtp)*e3t_n*tn + rbcp*e3t_b*tb )   ln_dynhpg_imp = T
+      !!                  /( rbcp*e3t_a    + (2-rbcp)*e3t_n    + rbcp*e3t_b    )   
+      !!             ztm = 0                                                       otherwise
       !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
       !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )
@@ -286,17 +284,14 @@
 
       !                                              ! ----------------------- !
-      IF( ln_dynhpg_imp ) THEN                       ! semi-implicite hpg case !
+      IF( ln_dynhpg_imp ) THEN                       ! semi-implicite hpg case !  (mean tracer computation)
          !                                           ! ----------------------- !
          !
          IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step
-            DO jk = 1, jpkm1                              ! (only swap)
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     tn(ji,jj,jk) = ta(ji,jj,jk)                    ! tn <-- ta
-                     sn(ji,jj,jk) = sa(ji,jj,jk)
-                  END DO
-               END DO
-            END DO
-         ELSE
+            DO jk = 1, jpkm1                                   ! (only swap)
+               tn(:,:,jk) = ta(:,:,jk)                         ! tn <-- ta
+               sn(:,:,jk) = sa(:,:,jk)
+            END DO
+            !                                             ! general case (Leapfrog + Asselin filter)
+         ELSE                                             ! apply filter on thickness weighted tracer and swap
             DO jk = 1, jpkm1
                DO jj = 1, jpj
@@ -311,7 +306,7 @@
                      !
                      !                                         ! Asselin filter on thickness and tracer content
-                     ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b )
-                     ztc_f  = atfp * ( ztca   - 2.* ztcn   + ztcb   ) 
-                     zsc_f  = atfp * ( zsca   - 2.* zscn   + zscb   ) 
+                     ze3t_f = atfp * ( ze3t_a + ze3t_b - 2.* ze3t_n )
+                     ztc_f  = atfp * ( ztca   + ztcb   - 2.* ztcn   ) 
+                     zsc_f  = atfp * ( zsca   + zscb   - 2.* zscn   ) 
                      !
                      !                                         ! filtered tracer including the correction 
@@ -320,9 +315,9 @@
                      zsf   = ze3fr * ( zscn   + zsc_f  )
                      !                                         ! mean thickness and tracer
-                     ze3mr = 1.e0  / ( ze3t_a + 2.* ze3t_n + ze3t_b )
-                     ztm   = ze3mr * ( ztca   + 2.* ztcn   + ztcb   )
-                     zsm   = ze3mr * ( zsca   + 2.* zscn   + zscb   )
+                     ze3mr = 1.e0  / (  rbcp * ( ze3t_a + ze3t_b ) + r1_2bcp * ze3t_n )
+                     ztm   = ze3mr * (  rbcp * ( ztca   + ztcb   ) + r1_2bcp * ztcn   )
+                     zsm   = ze3mr * (  rbcp * ( zsca   + zscb   ) + r1_2bcp * zscn   )
 !!gm mean e3t have to be saved and used in dynhpg  or it can be recomputed in dynhpg !!
-!!gm                 e3t_m(ji,jj,jk) = 0.25 / ze3mr
+!!gm                 e3t_m(ji,jj,jk) = 1.e0 / ze3mr
                      !                                         ! swap of arrays
                      tb(ji,jj,jk) = ztf                             ! tb <-- tn + filter
@@ -337,15 +332,11 @@
          ENDIF
          !                                           ! ----------------------- !
-      ELSE                                           !    explicit hpg case    !
+      ELSE                                           !    explicit hpg case    !  (NO mean tracer computation)
          !                                           ! ----------------------- !
          !
          IF( neuler == 0 .AND. kt == nit000 ) THEN        ! case of Euler time-stepping at first time-step
-            DO jk = 1, jpkm1                              ! No filter nor thickness weighting computation required
-               DO jj = 1, jpj                             ! ONLY swap
-                  DO ji = 1, jpi
-                     tn(ji,jj,jk) = ta(ji,jj,jk)                                 ! tn <-- ta
-                     sn(ji,jj,jk) = sa(ji,jj,jk)
-                  END DO
-               END DO
+            DO jk = 1, jpkm1                                   ! (only swap)
+               tn(:,:,jk) = ta(:,:,jk)                         ! tn <-- ta
+               sn(:,:,jk) = sa(:,:,jk)
             END DO
             !                                             ! general case (Leapfrog + Asselin filter)
