Index: branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/TRA/tranxt.F90
===================================================================
--- branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 2120)
+++ branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 2148)
@@ -15,14 +15,16 @@
    !!            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-05  (C. Ethe, G. Madec)  merge TRC-TRA
+   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  semi-implicit hpg with asselin filter + modified LF-RA
+   !!             -   !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
-   !!   tra_nxt       : time stepping on temperature and salinity
-   !!   tra_nxt_fix   : time stepping on temperature and salinity : fixed    volume case
-   !!   tra_nxt_vvl   : time stepping on temperature and salinity : variable volume case
+   !!   tra_nxt       : time stepping on tracers
+   !!   tra_nxt_fix   : time stepping on tracers : fixed    volume case
+   !!   tra_nxt_vvl   : time stepping on tracers : variable volume case
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers variables
    USE dom_oce         ! ocean space and time domain variables 
+   USE sbc_oce         ! surface boundary condition: ocean
    USE zdf_oce         ! ???
    USE domvvl          ! variable volume
@@ -38,4 +40,5 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE prtctl          ! Print control
+   USE traqsr          ! penetrative solar radiation (needed for nksr)
    USE traswp          ! swap array
    USE agrif_opa_update
@@ -49,4 +52,5 @@
    PUBLIC   tra_nxt_vvl   ! to be used in trcnxt
 
+   REAL(wp)                 ::   rbcp            ! Brown & Campana parameters for semi-implicit hpg
    REAL(wp), DIMENSION(jpk) ::   r2dt   ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler)
 
@@ -96,4 +100,6 @@
          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) * ( 1. - atfp)       ! Brown & Campana parameter for semi-implicit hpg
       ENDIF
 
@@ -107,5 +113,5 @@
 #endif
 
-#if defined key_obc
+#if defined key_obc 
       IF( lk_obc )   CALL obc_tra( kt )  ! OBC open boundaries
 #endif
@@ -114,5 +120,5 @@
 #endif
 #if defined key_agrif
-      CALL Agrif_tra                     ! AGRIF zoom boundaries
+      CALL Agrif_tra                   ! AGRIF zoom boundaries
 #endif
 
@@ -133,6 +139,6 @@
 
       ! Leap-Frog + Asselin filter time stepping
-      IF( lk_vvl )   THEN   ;   CALL tra_nxt_vvl( kt, tsb, tsn, tsa, jpts )  ! variable volume level (vvl)     
-      ELSE                  ;   CALL tra_nxt_fix( kt, tsb, tsn, tsa, jpts )  ! fixed    volume level 
+      IF( lk_vvl )   THEN   ;   CALL tra_nxt_vvl( kt, 'TRA', tsb, tsn, tsa, jpts )  ! variable volume level (vvl)     
+      ELSE                  ;   CALL tra_nxt_fix( kt,        tsb, tsn, tsa, jpts )  ! fixed    volume level 
       ENDIF
 
@@ -176,8 +182,9 @@
       !!              - 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
+      !!             ztm = tn + rbcp * [ta -2 tn + tb ]       ln_dynhpg_imp = T
+      !!             ztm = 0                                   otherwise
+      !!                   with rbcp=1/4 * (1-atfp^4) / (1-atfp)
       !!             tb  = tn + atfp*[ tb - 2 tn + ta ]
-      !!             tn  = ta 
+      !!             tn  = ta  
       !!             ta  = ztm       (NB: reset to 0 after eos_bn2 call)
       !!
@@ -185,12 +192,12 @@
       !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
       !!----------------------------------------------------------------------
-      INTEGER , INTENT(in   )                               ::  kt            ! ocean time-step index
-      INTEGER , INTENT(in   )                               ::  kjpt            ! number of tracers
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptb  ! before tracer fields
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptn  ! now tracer fields
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  pta         ! tracer trend
+      INTEGER , INTENT(in   )                               ::  kt       ! ocean time-step index
+      INTEGER , INTENT(in   )                               ::  kjpt     ! number of tracers
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptb      ! before tracer fields
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptn      ! now tracer fields
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  pta      ! tracer trend
       !!
       INTEGER  :: ji, jj, jk, jn   ! dummy loop indices
-      REAL(wp) :: ztm, ztf     ! temporary scalars
+      REAL(wp) :: ztd, ztm         ! temporary scalars
       !!----------------------------------------------------------------------
 
@@ -205,5 +212,5 @@
          !                                             ! (only swap)
          DO jn = 1, kjpt
-            DO jk = 1, jpkm1                              
+            DO jk = 1, jpkm1
                ptn(:,:,jk,jn) = pta(:,:,jk,jn)     ! ptb <-- ptn
             END DO
@@ -219,9 +226,10 @@
                   DO jj = 1, jpj
                      DO ji = 1, jpi
-                        ztm = 0.25 * ( pta(ji,jj,jk,jn) + 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn) )  ! mean pt
-                        ztf = atfp * ( pta(ji,jj,jk,jn) - 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn) )  ! Asselin filter on pt 
-                        ptb(ji,jj,jk,jn) = ptn(ji,jj,jk,jn) + ztf                                      ! ptb <-- filtered ptn 
-                        ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                                            ! ptn <-- pta
-                        pta(ji,jj,jk,jn) = ztm                                                           ! pta <-- mean pt
+                        ztd = pta(ji,jj,jk,jn) - 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn)    !  time laplacian on tracers
+                        ztm = ptn(ji,jj,jk,jn) + rbcp * ztd                                 !  used for both Asselin and Brown & Campana filters
+                        !
+                        ptb(ji,jj,jk,jn) = ptn(ji,jj,jk,jn) + atfp * ztd                    ! ptb <-- filtered ptn 
+                        ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                                 ! ptn <-- pta
+                        pta(ji,jj,jk,jn) = ztm                                              ! pta <-- Brown & Campana average
                      END DO
                   END DO
@@ -235,12 +243,12 @@
                   DO jj = 1, jpj
                      DO ji = 1, jpi
-                        ztf = atfp * ( pta(ji,jj,jk,jn) - 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn) )       ! Asselin filter on t 
-                        ptb(ji,jj,jk,jn) = ptn(ji,jj,jk,jn) + ztf                                     ! ptb <-- filtered ptn 
-                        ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                                           ! ptn <-- pta
+                        ztd = pta(ji,jj,jk,jn) - 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn)    ! time laplacian on tracers
+                        ! 
+                        ptb(ji,jj,jk,jn) = ptn(ji,jj,jk,jn) + atfp * ztd                    ! ptb <-- filtered ptn 
+                        ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                                 ! ptn <-- pta
                      END DO
                   END DO
                END DO
             END DO
-            !
          ENDIF
          !
@@ -250,5 +258,5 @@
 
 
-   SUBROUTINE tra_nxt_vvl( kt, ptb, ptn, pta, kjpt )
+   SUBROUTINE tra_nxt_vvl( kt, cdtype, ptb, ptn, pta, kjpt )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE tra_nxt_vvl  ***
@@ -263,7 +271,7 @@
       !!              - 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
+      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T
+      !!                  /( e3t_n    + rbcp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )   
+      !!             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    ] )
@@ -274,14 +282,17 @@
       !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
       !!----------------------------------------------------------------------
-      INTEGER , INTENT(in   )                               ::  kt            ! ocean time-step index
-      INTEGER , INTENT(in   )                               ::  kjpt            ! number of tracers
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptb  ! before tracer fields
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptn  ! now tracer fields
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  pta         ! tracer trend
+      INTEGER         , INTENT(in   )                               ::  kt       ! ocean time-step index
+      CHARACTER(len=3), INTENT(in   )                               ::  cdtype   ! =TRA or TRC (tracer indicator)
+      INTEGER         , INTENT(in   )                               ::  kjpt     ! number of tracers
+      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptb      ! before tracer fields
+      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptn      ! now tracer fields
+      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  pta      ! tracer trend
       !!     
-      INTEGER  ::   ji, jj, jk, jn             ! dummy loop indices
-      REAL(wp) ::   ztm , ztc_f , ztf , ztca, ztcn, ztcb   ! temporary scalar
-      REAL(wp) ::   ze3mr, ze3fr                           !    -         -
-      REAL(wp) ::   ze3t_b, ze3t_n, ze3t_a, ze3t_f         !    -         -
+      INTEGER  ::   ji, jj, jk, jn                 ! dummy loop indices
+      REAL(wp) ::   ztc_a , ztc_n , ztc_b          ! temporary scalar
+      REAL(wp) ::   ztc_f , ztc_d , ztc_m          !    -         -
+      REAL(wp) ::   ze3t_b, ze3t_n, ze3t_a         !    -         -
+      REAL(wp) ::   ze3t_f, ze3t_d, ze3t_m         !    -         -
+      REAL     ::   zfact1, zfact2                 !    -         -
       !!----------------------------------------------------------------------
 
@@ -293,18 +304,18 @@
       !
       !
-      IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step
-         !                                             ! (only swap)
+      IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step
+         !                                           ! (only swap)
          DO jn = 1, kjpt
-            DO jk = 1, jpkm1                              
-               ptn(:,:,jk,jn) = pta(:,:,jk,jn)         ! ptb <-- ptn
+            DO jk = 1, jpkm1
+               ptn(:,:,jk,jn) = pta(:,:,jk,jn)       ! ptb <-- ptn
             END DO
          END DO
          !                                              
-      ELSE                                              ! general case (Leapfrog + Asselin filter
+      ELSE                                           ! general case (Leapfrog + Asselin filter)
          !
-         !                                              ! ----------------------- !
-         IF( ln_dynhpg_imp ) THEN                       ! semi-implicite hpg case !
-            !                                           ! ----------------------- !
-            DO jn = 1, kjpt                               
+         !                                           ! ----------------------- !
+         IF( ln_dynhpg_imp ) THEN                    ! semi-implicite hpg case !
+            !                                        ! ----------------------- !
+            DO jn = 1, kjpt                          
                DO jk = 1, jpkm1
                   DO jj = 1, jpj
@@ -314,57 +325,89 @@
                         ze3t_a = fse3t_a(ji,jj,jk)
                         !                                         ! tracer content at Before, now and after
-                        ztcb = ptb(ji,jj,jk,jn) *  ze3t_b  
-                        ztcn = ptn(ji,jj,jk,jn) *  ze3t_n  
-                        ztca = pta(ji,jj,jk,jn) *  ze3t_a  
-                        !
-                        !                                         ! Asselin filter on thickness and tracer content
-                        ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b )
-                        ztc_f  = atfp * ( ztca   - 2.* ztcn   + ztcb   ) 
-                        !
-                        !                                         ! filtered tracer including the correction 
-                        ze3fr = 1.e0  / ( ze3t_n + ze3t_f )
-                        ztf   = ze3fr * ( ztcn   + ztc_f  )
-                        !                                         ! mean thickness and tracer
-                        ze3mr = 1.e0  / ( ze3t_a + 2.* ze3t_n + ze3t_b )
-                        ztm   = ze3mr * ( ztca   + 2.* ztcn   + ztcb   )
-                        !!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
+                        ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b  
+                        ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n  
+                        ztc_a  = pta(ji,jj,jk,jn) * ze3t_a  
+                        !                                         ! Time laplacian on tracer contents
+                        !                                         ! used for both Asselin and Brown & Campana filters
+                        ze3t_d =  ze3t_b - 2.* ze3t_n + ze3t_a 
+                        ztc_d  =  ztc_b  - 2.* ztc_n  + ztc_a  
+                        !                                         ! Asselin Filter on thicknesses and tracer contents
+                        ztc_f  = ztc_n  + atfp * ztc_d
+                        ztc_m  = ztc_n  + rbcp * ztc_d
+                        !                                         
+                        ze3t_f = 1.0 / ( ze3t_n + atfp * ze3t_d )
+                        ze3t_m = 1.0 / ( ze3t_n + rbcp * ze3t_d )
                         !                                         ! swap of arrays
-                        ptb(ji,jj,jk,jn) = ztf                             ! ptb <-- ptn + filter
-                        ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)              ! ptn <-- pta
-                        pta(ji,jj,jk,jn) = ztm                             ! pta <-- mean t
+                        ptb(ji,jj,jk,jn) = ztc_f * ze3t_f              ! ptb <-- ptn filtered
+                        pta(ji,jj,jk,jn) = ztc_m * ze3t_m              ! pta <-- Brown & Campana average
+                        ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)            ! ptn <-- pta
                      END DO
                   END DO
                END DO
             END DO
-            !                                        ! ----------------------- !
-         ELSE                                        !    explicit hpg case    !
-            !                                        ! ----------------------- !
-            DO jn = 1, kjpt      
-               DO jk = 1, jpkm1
-                  DO jj = 1, jpj
-                     DO ji = 1, jpi
-                        ze3t_b = fse3t_b(ji,jj,jk)
-                        ze3t_n = fse3t_n(ji,jj,jk)
-                        ze3t_a = fse3t_a(ji,jj,jk)
-                        !                                         ! tracer content at Before, now and after
-                        ztcb = ptb(ji,jj,jk,jn) *  ze3t_b  
-                        ztcn = ptn(ji,jj,jk,jn) *  ze3t_n   
-                        ztca = pta(ji,jj,jk,jn) *  ze3t_a  
-                        !
-                        !                                         ! Asselin filter on thickness and tracer content
-                        ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b )
-                        ztc_f  = atfp * ( ztca   - 2.* ztcn   + ztcb   ) 
-                        !
-                        !                                         ! filtered tracer including the correction 
-                        ze3fr = 1.e0 / ( ze3t_n + ze3t_f )
-                        ztf   =  ( ztcn  + ztc_f ) * ze3fr
-                        !                                         ! swap of arrays
-                        ptb(ji,jj,jk,jn) = ztf                    ! tb <-- tn filtered
-                        ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)       ! tn <-- ta
+            !                                           ! ----------------------- !
+         ELSE                                           !    explicit hpg case    !
+            !                                           ! ----------------------- !
+            IF( cdtype == 'TRA' ) THEN
+               !              
+               DO jn = 1, kjpt      
+                  DO jk = 1, jpkm1
+                     zfact1 = atfp * rdttra(jk)
+                     zfact2 = zfact1 / rau0
+                     DO jj = 1, jpj
+                        DO ji = 1, jpi
+                           ze3t_b = fse3t_b(ji,jj,jk)
+                           ze3t_n = fse3t_n(ji,jj,jk)
+                           ze3t_a = fse3t_a(ji,jj,jk)
+                           !                                         ! tracer content at Before, now and after
+                           ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b
+                           ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n
+                           ztc_a  = pta(ji,jj,jk,jn) * ze3t_a
+                           !
+                           ze3t_f = ze3t_n + atfp * ( ze3t_b - 2. * ze3t_n + ze3t_a )
+                           ztc_f  = ztc_n  + atfp * ( ztc_a  - 2. * ztc_n  + ztc_b )
+
+                           IF( jk == 1 ) THEN
+                               ze3t_f = ze3t_f - zfact2 * ( emp_b(ji,jj) - emp(ji,jj) )
+                               IF( jn == jp_tem )   ztc_f  = ztc_f  - zfact1 * ( sbc_hc_n(ji,jj) - sbc_hc_b(ji,jj) )
+                               IF( jn == jp_sal )   ztc_f  = ztc_f  - zfact1 * ( sbc_sc_n(ji,jj) - sbc_sc_b(ji,jj) )
+                           ENDIF
+                           IF( jn == jp_tem .AND. ln_traqsr .AND. jk <= nksr )  &
+                           &                        ztc_f  = ztc_f  - zfact1 * ( qsr_hc_n(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 
+
+                           ze3t_f = 1.e0 / ze3t_f
+                           ptb(ji,jj,jk,jn) = ztc_f * ze3t_f                           ! tb <-- tn filtered
+                           ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                             ! tn <-- ta
+                        END DO
                      END DO
                   END DO
                END DO
-            END DO
+               !
+            ELSE IF( cdtype == 'TRC' ) THEN
+               !              
+               DO jn = 1, kjpt      
+                  DO jk = 1, jpkm1
+                     DO jj = 1, jpj
+                        DO ji = 1, jpi
+                           ze3t_b = fse3t_b(ji,jj,jk)
+                           ze3t_n = fse3t_n(ji,jj,jk)
+                           ze3t_a = fse3t_a(ji,jj,jk)
+                           !                                         ! tracer content at Before, now and after
+                           ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b
+                           ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n
+                           ztc_a  = pta(ji,jj,jk,jn) * ze3t_a
+                           !
+                           ze3t_f = ze3t_n + atfp * ( ze3t_b - 2. * ze3t_n + ze3t_a )
+                           ztc_f  = ztc_n  + atfp * ( ztc_a  - 2. * ztc_n  + ztc_b  )
+
+                           ze3t_f = 1.e0 / ze3t_f
+                           ptb(ji,jj,jk,jn) = ztc_f * ze3t_f                           ! tb <-- tn filtered
+                           ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                             ! tn <-- ta
+                        END DO
+                     END DO
+                  END DO
+               END DO
+               !
+            END IF
             !
          ENDIF
