Index: /NEMO/branches/UKMO/NEMO_4.0.4_CO9_OBS_ASM/src/OCE/ASM/asmbkg.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.4_CO9_OBS_ASM/src/OCE/ASM/asmbkg.F90	(revision 15182)
+++ /NEMO/branches/UKMO/NEMO_4.0.4_CO9_OBS_ASM/src/OCE/ASM/asmbkg.F90	(revision 15183)
@@ -39,4 +39,5 @@
    USE ice
 #endif
+   USE asminc, ONLY: ln_avgbkg
 
    IMPLICIT NONE
@@ -44,4 +45,12 @@
    
    PUBLIC   asm_bkg_wri   !: Write out the background state
+
+  !! * variables for calculating time means
+   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   tn_tavg  , sn_tavg  
+   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   un_tavg  , vn_tavg
+   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   avt_tavg
+   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   en_tavg
+   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:)   ::   sshn_tavg
+   REAL(wp),SAVE :: numtimes_tavg     ! No of times to average over
 
    !!----------------------------------------------------------------------
@@ -71,9 +80,73 @@
       INTEGER :: inum          ! File unit number
       REAL(wp) :: zdate        ! Date
+      INTEGER :: ierror
       !!-----------------------------------------------------------------------
 
-      !                                !-------------------------------------------
-      IF( kt == nitbkg_r ) THEN        ! Write out background at time step nitbkg_r
-         !                             !-----------------------------------========
+      ! If creating an averaged assim bkg, initialise on first timestep
+      IF ( ln_avgbkg .AND. kt == ( nn_it000 - 1) ) THEN
+         ! Allocate memory 
+         ALLOCATE( tn_tavg(jpi,jpj,jpk), STAT=ierror )
+         IF( ierror > 0 ) THEN
+            CALL ctl_stop( 'asm_bkg_wri: unable to allocate tn_tavg' )   ;   RETURN
+         ENDIF
+         tn_tavg(:,:,:)=0
+         ALLOCATE( sn_tavg(jpi,jpj,jpk), STAT=ierror )
+         IF( ierror > 0 ) THEN
+            CALL ctl_stop( 'asm_bkg_wri: unable to allocate sn_tavg' )   ;   RETURN
+         ENDIF
+         sn_tavg(:,:,:)=0
+         ALLOCATE( un_tavg(jpi,jpj,jpk), STAT=ierror )
+         IF( ierror > 0 ) THEN
+            CALL ctl_stop( 'asm_bkg_wri: unable to allocate un_tavg' )   ;   RETURN
+         ENDIF
+         un_tavg(:,:,:)=0
+         ALLOCATE( vn_tavg(jpi,jpj,jpk), STAT=ierror )
+         IF( ierror > 0 ) THEN
+            CALL ctl_stop( 'asm_bkg_wri: unable to allocate vn_tavg' )   ;   RETURN
+         ENDIF
+         vn_tavg(:,:,:)=0
+         ALLOCATE( sshn_tavg(jpi,jpj), STAT=ierror )
+         IF( ierror > 0 ) THEN
+            CALL ctl_stop( 'asm_bkg_wri: unable to allocate sshn_tavg' )   ;   RETURN
+         ENDIF
+         sshn_tavg(:,:)=0
+         IF( ln_zdftke ) THEN
+            ALLOCATE( en_tavg(jpi,jpj,jpk), STAT=ierror )
+            IF( ierror > 0 ) THEN
+               CALL ctl_stop( 'asm_bkg_wri: unable to allocate en_tavg' )   ;   RETURN
+            ENDIF
+            en_tavg(:,:,:)=0
+         ENDIF
+         ALLOCATE( avt_tavg(jpi,jpj,jpk), STAT=ierror )
+         IF( ierror > 0 ) THEN
+            CALL ctl_stop( 'asm_bkg_wri: unable to allocate avt_tavg' )   ;   RETURN
+         ENDIF
+         avt_tavg(:,:,:)=0
+         
+         numtimes_tavg = REAL ( nitavgbkg_r -  nn_it000 + 1 )
+      ENDIF
+      
+      ! If creating an averaged assim bkg, sum the contribution every timestep
+      IF ( ln_avgbkg ) THEN 
+         IF (lwp) THEN
+              WRITE(numout,*) 'asm_bkg_wri : Summing assim bkg fields at timestep ',kt
+              WRITE(numout,*) '~~~~~~~~~~~~ '
+         ENDIF
+
+         tn_tavg(:,:,:)  = tn_tavg(:,:,:)  + tsn(:,:,:,jp_tem) / numtimes_tavg
+         sn_tavg(:,:,:)  = sn_tavg(:,:,:)  + tsn(:,:,:,jp_sal) / numtimes_tavg
+         sshn_tavg(:,:)  = sshn_tavg(:,:)  + sshn (:,:)        / numtimes_tavg
+         un_tavg(:,:,:)  = un_tavg(:,:,:)  + un(:,:,:)         / numtimes_tavg
+         vn_tavg(:,:,:)  = vn_tavg(:,:,:)  + vn(:,:,:)         / numtimes_tavg
+         avt_tavg(:,:,:) = avt_tavg(:,:,:) + avt(:,:,:)        / numtimes_tavg
+         IF( ln_zdftke ) THEN
+            en_tavg(:,:,:)  = en_tavg(:,:,:)  + en(:,:,:)         / numtimes_tavg
+         ENDIF
+      ENDIF
+     
+
+      ! Write out background at time step nitbkg_r or nitavgbkg_r
+      IF ( ( .NOT. ln_avgbkg .AND. (kt == nitbkg_r) ) .OR. &
+      &          ( ln_avgbkg .AND. (kt == nitavgbkg_r) ) ) THEN
          !
          WRITE(cl_asmbkg, FMT='(A,".nc")' ) TRIM( c_asmbkg )
@@ -87,25 +160,45 @@
             CALL iom_open( c_asmbkg, inum, ldwrt = .TRUE. )
             !
-            IF( nitbkg_r == nit000 - 1 ) THEN      ! Treat special case when nitbkg = 0
-               zdate = REAL( ndastp )
-               IF( ln_zdftke ) THEN                   ! read turbulent kinetic energy ( en )
-                  IF(lwp) WRITE(numout,*) ' Reading TKE (en) from restart...'
-                  CALL tke_rst( nit000, 'READ' )
+            !                                      ! Write the information
+            IF ( ln_avgbkg ) THEN
+               IF( nitavgbkg_r == nit000 - 1 ) THEN      ! Treat special case when nitavgbkg = 0
+                  zdate = REAL( ndastp )
+                  IF( ln_zdftke ) THEN
+                     IF(lwp) WRITE(numout,*) ' Reading TKE (en) from restart...'
+                     CALL tke_rst( nit000, 'READ' )               ! lk_zdftke=T :   Read turbulent kinetic energy ( en )
+                  ENDIF
+               ELSE
+                  zdate = REAL( ndastp )
                ENDIF
+               CALL iom_rstput( kt, nitavgbkg_r, inum, 'rdastp' , zdate   )
+               CALL iom_rstput( kt, nitavgbkg_r, inum, 'un'     , un_tavg )
+               CALL iom_rstput( kt, nitavgbkg_r, inum, 'vn'     , vn_tavg )
+               CALL iom_rstput( kt, nitavgbkg_r, inum, 'tn'     , tn_tavg )
+               CALL iom_rstput( kt, nitavgbkg_r, inum, 'sn'     , sn_tavg )
+               CALL iom_rstput( kt, nitavgbkg_r, inum, 'sshn'   , sshn_tavg)
+               IF( ln_zdftke ) CALL iom_rstput( kt, nitavgbkg_r, inum, 'en'     , en_tavg )
+               CALL iom_rstput( kt, nitavgbkg_r, inum, 'avt'    , avt_tavg)
+               !
             ELSE
-               zdate = REAL( ndastp )
-            ENDIF
-            !
-            !                                      ! Write the information
-            CALL iom_rstput( kt, nitbkg_r, inum, 'rdastp' , zdate             )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'un'     , un                )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'vn'     , vn                )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'tn'     , tsn(:,:,:,jp_tem) )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'sn'     , tsn(:,:,:,jp_sal) )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'sshn'   , sshn              )
-            IF( ln_zdftke )   CALL iom_rstput( kt, nitbkg_r, inum, 'en'     , en                )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'avt'    , avt               )
-            !
-            CALL iom_close( inum )
+               IF( nitbkg_r == nit000 - 1 ) THEN      ! Treat special case when nitbkg = 0
+                  zdate = REAL( ndastp )
+                  IF( ln_zdftke ) THEN
+                     IF(lwp) WRITE(numout,*) ' Reading TKE (en) from restart...'
+                     CALL tke_rst( nit000, 'READ' )               ! lk_zdftke=T :   Read turbulent kinetic energy ( en )
+                  ENDIF
+               ELSE
+                  zdate = REAL( ndastp )
+               ENDIF
+               CALL iom_rstput( kt, nitbkg_r, inum, 'rdastp' , zdate             )
+               CALL iom_rstput( kt, nitbkg_r, inum, 'un'     , un                )
+               CALL iom_rstput( kt, nitbkg_r, inum, 'vn'     , vn                )
+               CALL iom_rstput( kt, nitbkg_r, inum, 'tn'     , tsn(:,:,:,jp_tem) )
+               CALL iom_rstput( kt, nitbkg_r, inum, 'sn'     , tsn(:,:,:,jp_sal) )
+               CALL iom_rstput( kt, nitbkg_r, inum, 'sshn'   , sshn              )
+               IF( ln_zdftke )  CALL iom_rstput( kt, nitbkg_r, inum, 'en'     , en                )
+               CALL iom_rstput( kt, nitbkg_r, inum, 'avt'    , avt               )
+               !
+               CALL iom_close( inum )
+            ENDIF
          ENDIF
          !
Index: /NEMO/branches/UKMO/NEMO_4.0.4_CO9_OBS_ASM/src/OCE/ASM/asminc.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.4_CO9_OBS_ASM/src/OCE/ASM/asminc.F90	(revision 15182)
+++ /NEMO/branches/UKMO/NEMO_4.0.4_CO9_OBS_ASM/src/OCE/ASM/asminc.F90	(revision 15183)
@@ -31,5 +31,4 @@
    USE zpshde          ! Partial step : Horizontal Derivative
    USE asmpar          ! Parameters for the assmilation interface
-   USE asmbkg          ! 
    USE c1d             ! 1D initialization
    USE sbc_oce         ! Surface boundary condition variables.
@@ -59,4 +58,5 @@
 #endif
    LOGICAL, PUBLIC :: ln_bkgwri     !: No output of the background state fields
+   LOGICAL, PUBLIC :: ln_avgbkg     !: No output of the mean background state fields
    LOGICAL, PUBLIC :: ln_asmiau     !: No applying forcing with an assimilation increment
    LOGICAL, PUBLIC :: ln_asmdin     !: No direct initialization
@@ -82,4 +82,5 @@
    INTEGER , PUBLIC ::   nitiaustr   !: Time step of the start of the IAU interval 
    INTEGER , PUBLIC ::   nitiaufin   !: Time step of the end of the IAU interval
+   INTEGER , PUBLIC ::   nitavgbkg   !: Number of timesteps to average assim bkg [0,nitavgbkg]
    ! 
    INTEGER , PUBLIC ::   niaufn      !: Type of IAU weighing function: = 0   Constant weighting
@@ -122,4 +123,5 @@
       REAL(KIND=dp) :: ditiaustr_date  ! Date YYYYMMDD.HHMMSS of IAU interval start time step
       REAL(KIND=dp) :: ditiaufin_date  ! Date YYYYMMDD.HHMMSS of IAU interval final time step
+      REAL(KIND=dp) :: ditavgbkg_date  ! Date YYYYMMDD.HHMMSS of end of assim bkg averaging period
 
       REAL(wp) :: znorm        ! Normalization factor for IAU weights
@@ -132,9 +134,9 @@
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zhdiv   ! 2D workspace
       !!
-      NAMELIST/nam_asminc/ ln_bkgwri,                                      &
+      NAMELIST/nam_asminc/ ln_bkgwri, ln_avgbkg,                           &
          &                 ln_trainc, ln_dyninc, ln_sshinc,                &
          &                 ln_asmdin, ln_asmiau,                           &
          &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   &
-         &                 ln_salfix, salfixmin, nn_divdmp
+         &                 ln_salfix, salfixmin, nn_divdmp, nitavgbkg
       !!----------------------------------------------------------------------
 
@@ -142,6 +144,22 @@
       ! Read Namelist nam_asminc : assimilation increment interface
       !-----------------------------------------------------------------------
-      ln_seaiceinc   = .FALSE.
+      ! Set default values
+      ln_bkgwri = .FALSE.
+      ln_avgbkg = .FALSE.
+      ln_trainc = .FALSE.
+      ln_dyninc = .FALSE.
+      ln_sshinc = .FALSE.
+      ln_asmdin = .FALSE.
+      ln_asmiau = .TRUE.
+      ln_salfix = .FALSE.
+      ln_seaiceinc = .FALSE.
       ln_temnofreeze = .FALSE.
+      nitbkg    = 0
+      nitdin    = 0      
+      nitiaustr = 1
+      nitiaufin = 150
+      niaufn    = 0
+      salfixmin = -9999
+      nitavgbkg = 1
 
       REWIND( numnam_ref )              ! Namelist nam_asminc in reference namelist : Assimilation increment
@@ -160,4 +178,5 @@
          WRITE(numout,*) '   Namelist namasm : set assimilation increment parameters'
          WRITE(numout,*) '      Logical switch for writing out background state          ln_bkgwri = ', ln_bkgwri
+         WRITE(numout,*) '      Logical switch for writing mean background state         ln_avgbkg = ', ln_avgbkg
          WRITE(numout,*) '      Logical switch for applying tracer increments            ln_trainc = ', ln_trainc
          WRITE(numout,*) '      Logical switch for applying velocity increments          ln_dyninc = ', ln_dyninc
@@ -170,4 +189,5 @@
          WRITE(numout,*) '      Timestep of start of IAU interval in [0,nitend-nit000-1] nitiaustr = ', nitiaustr
          WRITE(numout,*) '      Timestep of end of IAU interval in [0,nitend-nit000-1]   nitiaufin = ', nitiaufin
+         WRITE(numout,*) '      Number of timesteps to average assim bkg [0,nitavgbkg]   nitavgbkg = ', nitavgbkg
          WRITE(numout,*) '      Type of IAU weighting function                           niaufn    = ', niaufn
          WRITE(numout,*) '      Logical switch for ensuring that the sa > salfixmin      ln_salfix = ', ln_salfix
@@ -179,4 +199,5 @@
       nitiaustr_r = nitiaustr + nit000 - 1            ! Start of IAU interval referenced to nit000
       nitiaufin_r = nitiaufin + nit000 - 1            ! End of IAU interval referenced to nit000
+      nitavgbkg_r = nitavgbkg + nit000 - 1            ! Averaging period referenced to nit000
 
       iiauper     = nitiaufin_r - nitiaustr_r + 1     ! IAU interval length
@@ -188,4 +209,5 @@
       CALL calc_date( nitiaustr_r, ditiaustr_date )   ! IAU start time referenced to ndate0
       CALL calc_date( nitiaufin_r, ditiaufin_date )   ! IAU end time referenced to ndate0
+      CALL calc_date( nitavgbkg_r, ditavgbkg_date )   ! End of assim bkg averaging period referenced to ndate0
 
       IF(lwp) THEN
@@ -199,4 +221,5 @@
          WRITE(numout,*) '       nitiaustr_r = ', nitiaustr_r
          WRITE(numout,*) '       nitiaufin_r = ', nitiaufin_r
+         WRITE(numout,*) '       nitavgbkg_r = ', nitavgbkg_r
          WRITE(numout,*)
          WRITE(numout,*) '   Dates referenced to current cycle:'
@@ -209,4 +232,5 @@
          WRITE(numout,*) '       ditiaustr_date = ', ditiaustr_date
          WRITE(numout,*) '       ditiaufin_date = ', ditiaufin_date
+         WRITE(numout,*) '       ditavgbkg_date = ', ditavgbkg_date
       ENDIF
 
@@ -249,4 +273,9 @@
          &                ' the cycle interval')
 
+      IF ( nitavgbkg_r > nitend ) &
+         & CALL ctl_stop( ' nitavgbkg_r :',  &
+         &                ' Assim bkg averaging period is outside', &
+         &                ' the cycle interval')
+      
       IF ( nstop > 0 ) RETURN       ! if there are any errors then go no further
 
@@ -494,5 +523,5 @@
       !
       IF( lk_asminc ) THEN                            !==  data assimilation  ==!
-         IF( ln_bkgwri )   CALL asm_bkg_wri( nit000 - 1 )      ! Output background fields
+!         IF( ln_bkgwri )   CALL asm_bkg_wri( nit000 - 1 )      ! Output background fields
          IF( ln_asmdin ) THEN                                  ! Direct initialization
             IF( ln_trainc )   CALL tra_asm_inc( nit000 - 1 )      ! Tracers
Index: /NEMO/branches/UKMO/NEMO_4.0.4_CO9_OBS_ASM/src/OCE/ASM/asmpar.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.4_CO9_OBS_ASM/src/OCE/ASM/asmpar.F90	(revision 15182)
+++ /NEMO/branches/UKMO/NEMO_4.0.4_CO9_OBS_ASM/src/OCE/ASM/asmpar.F90	(revision 15183)
@@ -21,4 +21,5 @@
    INTEGER, PUBLIC ::   nitiaufin_r   !: IAU final time step referenced to nit000
    INTEGER, PUBLIC ::   nittrjfrq     !: Frequency of trajectory output for 4D-VAR
+   INTEGER, PUBLIC ::   nitavgbkg_r   !: Averaging period for assim bkg referenced to nit000
 
    !!----------------------------------------------------------------------
Index: /NEMO/branches/UKMO/NEMO_4.0.4_CO9_OBS_ASM/src/OCE/nemogcm.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.4_CO9_OBS_ASM/src/OCE/nemogcm.F90	(revision 15182)
+++ /NEMO/branches/UKMO/NEMO_4.0.4_CO9_OBS_ASM/src/OCE/nemogcm.F90	(revision 15183)
@@ -166,4 +166,5 @@
       END DO
 #else
+      IF( lk_asminc .AND. ln_bkgwri ) CALL asm_bkg_wri( nit000 - 1 )    ! Output background fields
       !
 # if defined key_agrif
Index: /NEMO/branches/UKMO/NEMO_4.0.4_CO9_OBS_ASM/src/OCE/step.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.4_CO9_OBS_ASM/src/OCE/step.F90	(revision 15182)
+++ /NEMO/branches/UKMO/NEMO_4.0.4_CO9_OBS_ASM/src/OCE/step.F90	(revision 15183)
@@ -296,4 +296,5 @@
 #endif
 
+      IF( ln_bkgwri )    CALL asm_bkg_wri  ( kstp )      ! output background fields
       IF( ln_diaobs  )   CALL dia_obs      ( kstp )      ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
 
