Index: branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 2111)
+++ branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 2148)
@@ -4,6 +4,7 @@
    !! Surface module :  provide to the ocean its surface boundary condition
    !!======================================================================
-   !! History :  3.0   !  07-2006  (G. Madec)  Original code
-   !!             -    !  08-2008  (S. Masson, E. .... ) coupled interface
+   !! History :  3.0  !  2006-07  (G. Madec)  Original code
+   !!            3.1  !  2008-08  (S. Masson, E. Maisonnave, G. Madec) coupled interface
+   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
    !!----------------------------------------------------------------------
 
@@ -49,5 +50,5 @@
 #  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)
@@ -86,5 +87,5 @@
 !!gm  IF( lk_sbc_cpl       ) THEN   ;   ln_cpl      = .TRUE.   ;   ELSE   ;   ln_cpl      = .FALSE.   ;   ENDIF
 
-      IF ( Agrif_Root() ) THEN
+      IF( Agrif_Root() ) THEN
         IF( lk_lim2 )            nn_ice      = 2
         IF( lk_lim3 )            nn_ice      = 3
@@ -179,6 +180,8 @@
       !!                CAUTION : never mask the surface stress field (tke sbc)
       !!
-      !! ** Action  : - set the ocean surface boundary condition, i.e.  
-      !!                utau, vtau, qns, qsr, emp, emps, qrp, erp
+      !! ** Action  : - set the ocean surface boundary condition at before and now 
+      !!                time step, i.e.  
+      !!                utau_b, vtau_b, qns_b, qsr_b, emp_n, emps_b, qrp_b, erp_b
+      !!                utau  , vtau  , qns  , qsr  , emp  , emps  , qrp  , erp
       !!              - updte the ice fraction : fr_i
       !!----------------------------------------------------------------------
@@ -186,17 +189,28 @@
       !!---------------------------------------------------------------------
 
-      CALL iom_setkt( kt + nn_fsbc - 1 )         !  in sbc, iom_put is called every nn_fsbc time step
-      !
-      ! ocean to sbc mean sea surface variables (ss._m)
-      ! ---------------------------------------
-      CALL sbc_ssm( kt )                         ! sea surface mean currents (at U- and V-points), 
-      !                                          ! temperature and salinity (at T-point) over nf_sbc time-step
-      !                                          ! (i.e. sst_m, sss_m, ssu_m, ssv_m)
-
-      ! sbc formulation
-      ! ---------------
-         
-      SELECT CASE( nsbc )                        ! Compute ocean surface boundary condition
-      !                                          ! (i.e. utau,vtau, qns, qsr, emp, emps)
+      !                                            ! ---------------------------------------- !
+      IF( kt /= nit000 ) THEN                      !          Swap of forcing fields          !
+         !                                         ! ---------------------------------------- !
+         utau_b(:,:) = utau(:,:)                         ! Swap the ocean forcing fields
+         vtau_b(:,:) = vtau(:,:)                         ! (except at nit000 where before fields
+         qns_b (:,:) = qns (:,:)                         !  are set at the end of the routine)
+         ! The 3D heat content due to qsr forcing is treated in traqsr
+         ! qsr_b (:,:) = qsr (:,:)
+         emp_b (:,:) = emp (:,:)
+         emps_b(:,:) = emps(:,:)
+      ENDIF
+      !                                            ! ---------------------------------------- !
+      !                                            !        forcing field computation         !
+      !                                            ! ---------------------------------------- !
+
+      CALL iom_setkt( kt + nn_fsbc - 1 )                 ! in sbc, iom_put is called every nn_fsbc time step
+      !
+      CALL sbc_ssm( kt )                                 ! ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m)
+      !                                                  ! averaged over nf_sbc time-step
+
+                                                   !==  sbc formulation  ==!
+                                                            
+      SELECT CASE( nsbc )                                ! Compute ocean surface boundary condition
+      !                                                  ! (i.e. utau,vtau, qns, qsr, emp, emps)
       CASE(  0 )   ;   CALL sbc_gyre    ( kt )                    ! analytical formulation : GYRE configuration
       CASE(  1 )   ;   CALL sbc_ana     ( kt )                    ! analytical formulation : uniform sbc
@@ -214,6 +228,5 @@
       END SELECT
 
-      ! Misc. Options
-      ! -------------
+      !                                            !==  Misc. Options  ==!
 
 !!gm  IF( ln_dm2dc       )   CALL sbc_dcy( kt )                 ! Daily mean qsr distributed over the Diurnal Cycle
@@ -236,11 +249,52 @@
       !                                                         ! (update freshwater fluxes)
       !
+      IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    !
+         !                                             ! ---------------------------------------- !
+         IF( ln_rstart .AND.    &                               !* Restart: read in restart file
+            & iom_varid( numror, 'utau_b', ldstop = .FALSE. ) > 0 ) THEN 
+            IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields red in the restart file'
+            CALL iom_get( numror, jpdom_autoglo, 'utau_b', utau_b )   ! before i-stress  (U-point)
+            CALL iom_get( numror, jpdom_autoglo, 'vtau_b', vtau_b )   ! before j-stress  (V-point)
+            CALL iom_get( numror, jpdom_autoglo, 'qns_b' , qns_b  )   ! before non solar heat flux (T-point)
+            ! The 3D heat content due to qsr forcing is treated in traqsr
+            ! CALL iom_get( numror, jpdom_autoglo, 'qsr_b' , qsr_b  )   ! before     solar heat flux (T-point)
+            CALL iom_get( numror, jpdom_autoglo, 'emp_b' , emp_b  )   ! before     freshwater flux (T-point)
+            CALL iom_get( numror, jpdom_autoglo, 'emps_b', emps_b )   ! before C/D freshwater flux (T-point)
+         ELSE                                                   !* no restart: set from nit000 values
+            IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields set to nit000'
+            utau_b(:,:) = utau(:,:) 
+            vtau_b(:,:) = vtau(:,:)
+            qns_b (:,:) = qns (:,:)
+            ! qsr_b (:,:) = qsr (:,:)
+            emp_b (:,:) = emp (:,:)
+            emps_b(:,:) = emps(:,:)
+         ENDIF
+      ENDIF
+      !                                                ! ---------------------------------------- !
+      IF( lrst_oce ) THEN                              !      Write in the ocean restart file     !
+         !                                             ! ---------------------------------------- !
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'sbc : ocean surface forcing fields written in ocean restart file ',   &
+            &                    'at it= ', kt,' date= ', ndastp
+         IF(lwp) WRITE(numout,*) '~~~~'
+         CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , utau )
+         CALL iom_rstput( kt, nitrst, numrow, 'vtau_b' , vtau )
+         CALL iom_rstput( kt, nitrst, numrow, 'qns_b'  , qns  )
+         ! The 3D heat content due to qsr forcing is treated in traqsr
+         ! CALL iom_rstput( kt, nitrst, numrow, 'qsr_b'  , qsr  )
+         CALL iom_rstput( kt, nitrst, numrow, 'emp_b'  , emp  )
+         CALL iom_rstput( kt, nitrst, numrow, 'emps_b' , emps )
+      ENDIF
+
+      !                                                ! ---------------------------------------- !
+      !                                                !        Outputs and control print         !
+      !                                                ! ---------------------------------------- !
       IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
-         CALL iom_put( "emp-rnf"  , (emp-rnf)  )                ! upward water flux 
-         CALL iom_put( "emps-rnf" , (emps-rnf) )                ! c/d water flux 
-         CALL iom_put( "qns+qsr"  , qns + qsr  )                ! total heat flux   (caution if ln_dm2dc=true, to be 
-         CALL iom_put( "qns"      , qns        )                ! solar heat flux    moved after the call to iom_setkt)
-         CALL iom_put( "qsr"      ,       qsr  )                ! solar heat flux    moved after the call to iom_setkt)
-         IF(  nn_ice > 0 )   CALL iom_put( "ice_cover", fr_i )  ! ice fraction 
+         CALL iom_put( "emp-rnf" , emp  - rnf )                   ! upward water flux
+         CALL iom_put( "emps-rnf", emps - rnf )                   ! c/d water flux
+         CALL iom_put( "qns+qsr" , qns  + qsr )                   ! total heat flux   (caution if ln_dm2dc=true, to be 
+         CALL iom_put( "qns"     , qns        )                   ! solar heat flux    moved after the call to iom_setkt)
+         CALL iom_put( "qsr"     ,       qsr  )                   ! solar heat flux    moved after the call to iom_setkt)
+         IF( nn_ice > 0 )   CALL iom_put( "ice_cover", fr_i )   ! ice fraction 
       ENDIF
       !
@@ -253,14 +307,14 @@
       !
       IF(ln_ctl) THEN         ! print mean trends (used for debugging)
-         CALL prt_ctl(tab2d_1=fr_i   , clinfo1=' fr_i - : ', mask1=tmask, ovlap=1 )
-         CALL prt_ctl(tab2d_1=(emp-rnf) , clinfo1=' emp-rnf  - : ', mask1=tmask, ovlap=1 ) 
-         CALL prt_ctl(tab2d_1=(emps-rnf), clinfo1=' emps-rnf - : ', mask1=tmask, ovlap=1 ) 
-         CALL prt_ctl(tab2d_1=qns    , clinfo1=' qns  - : ', mask1=tmask, ovlap=1 )
-         CALL prt_ctl(tab2d_1=qsr    , clinfo1=' qsr  - : ', mask1=tmask, ovlap=1 )
-         CALL prt_ctl(tab3d_1=tmask  , clinfo1=' tmask  : ', mask1=tmask, ovlap=1, kdim=jpk )
-         CALL prt_ctl(tab3d_1=tn     , clinfo1=' sst  - : ', mask1=tmask, ovlap=1, kdim=1   )
-         CALL prt_ctl(tab3d_1=sn     , clinfo1=' sss  - : ', mask1=tmask, ovlap=1, kdim=1   )
-         CALL prt_ctl(tab2d_1=utau   , clinfo1=' utau - : ', mask1=umask,                      &
-            &         tab2d_2=vtau   , clinfo2=' vtau - : ', mask2=vmask, ovlap=1 )
+         CALL prt_ctl(tab2d_1=fr_i      , clinfo1=' fr_i     - : ', mask1=tmask, ovlap=1 )
+         CALL prt_ctl(tab2d_1=(emp-rnf) , clinfo1=' emp-rnf  - : ', mask1=tmask, ovlap=1 )
+         CALL prt_ctl(tab2d_1=(emps-rnf), clinfo1=' emps-rnf - : ', mask1=tmask, ovlap=1 )
+         CALL prt_ctl(tab2d_1=qns       , clinfo1=' qns      - : ', mask1=tmask, ovlap=1 )
+         CALL prt_ctl(tab2d_1=qsr       , clinfo1=' qsr      - : ', mask1=tmask, ovlap=1 )
+         CALL prt_ctl(tab3d_1=tmask     , clinfo1=' tmask    - : ', mask1=tmask, ovlap=1, kdim=jpk )
+         CALL prt_ctl(tab3d_1=tn        , clinfo1=' sst      - : ', mask1=tmask, ovlap=1, kdim=1   )
+         CALL prt_ctl(tab3d_1=sn        , clinfo1=' sss      - : ', mask1=tmask, ovlap=1, kdim=1   )
+         CALL prt_ctl(tab2d_1=utau      , clinfo1=' utau     - : ', mask1=umask,                      &
+            &         tab2d_2=vtau      , clinfo2=' vtau     - : ', mask2=vmask, ovlap=1 )
       ENDIF
       !
