Index: /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 8399)
+++ /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 8400)
@@ -1285,4 +1285,6 @@
     ln_asmdin = .false.    !  Logical switch for Direct Initialization (DI)
     ln_asmiau = .false.    !  Logical switch for Incremental Analysis Updating (IAU)
+    ln_seaiceinc = .false. !  Logical switch for applying sea ice increments
+    ln_temnofreeze = .false. !  Logical to not add increments if temperature would fall below freezing
     nitbkg    = 0          !  Timestep of background in [0,nitend-nit000-1]
     nitdin    = 0          !  Timestep of background for DI in [0,nitend-nit000-1]
@@ -1321,2 +1323,35 @@
    rn_htrmax         =  200.0   ! max. depth of transition range
 /
+!-----------------------------------------------------------------------
+&nambias   ! Bias pressure correctiom
+!-----------------------------------------------------------------------
+   ln_bias        = .false.
+   ln_bias_asm    = .false.
+   ln_bias_rlx    = .false.
+   ln_bias_ofl    = .false.
+   ln_bias_ts_app = .false.
+   ln_bias_pc_app = .false.        
+   fb_t_asm       = 0.0
+   fb_t_rlx       = 0.0
+   fb_t_ofl       = 1.0
+   fb_p_asm       = 1.0
+   fb_p_rlx       = 1.0
+   fb_p_ofl       = 0.0
+   eft_rlx        = 365.0
+   eft_asm        = 365.0
+   t_rlx_upd      = 0.1
+   t_asm_upd      = 0.1
+   nn_lat_ramp    = 0          
+   bias_time_unit_asm = 86400.0
+   bias_time_unit_rlx = 1.0
+   bias_time_unit_ofl = 1.0 
+   cn_bias_tot    = "bias_tot.nc" 
+   cn_bias_asm    = "bias_asm.nc"
+   cn_dir         = './'  
+   ln_bsyncro     = .FALSE. 
+   fctamp         = 1.
+   rn_maxlat_bias = 23.0      
+   rn_minlat_bias = 10.0
+   nn_bias_itwrt  = 15
+   ln_itdecay     = .FALSE.
+/
Index: /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90	(revision 8399)
+++ /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90	(revision 8400)
@@ -119,5 +119,6 @@
             CALL iom_rstput( kt, nitbkg_r, inum, 'en'     , en                )
 #endif
-            CALL iom_rstput( kt, nitbkg_r, inum, 'gcx'    , gcx               )
+!            CALL iom_rstput( kt, nitbkg_r, inum, 'gcx'    , gcx               )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'avt'    , avt               )
             !
             CALL iom_close( inum )
@@ -153,4 +154,5 @@
             CALL iom_rstput( kt, nitdin_r, inum, 'tn'     , tsn(:,:,:,jp_tem) )
             CALL iom_rstput( kt, nitdin_r, inum, 'sn'     , tsn(:,:,:,jp_sal) )
+            CALL iom_rstput( kt, nitdin_r, inum, 'avt'    , avt     )
             CALL iom_rstput( kt, nitdin_r, inum, 'sshn'   , sshn              )
 #if defined key_lim2 || defined key_lim3
Index: /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 8399)
+++ /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 8400)
@@ -39,4 +39,8 @@
    USE ice_2            ! LIM2
 #endif
+#if defined key_cice && defined key_asminc
+   USE sbc_ice, ONLY : & ! CICE Ice model variables
+   & ndaice_da, nfresh_da, nfsalt_da
+#endif
    USE sbc_oce          ! Surface boundary condition variables.
 
@@ -133,5 +137,6 @@
          &                 ln_asmdin, ln_asmiau,                           &
          &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   &
-         &                 ln_salfix, salfixmin, nn_divdmp
+         &                 ln_salfix, salfixmin, nn_divdmp,                &
+         &                 ln_seaiceinc, ln_temnofreeze
       !!----------------------------------------------------------------------
 
@@ -892,4 +897,8 @@
             ENDIF
 
+         ELSE
+#if defined key_asminc
+            ssh_iau(:,:) = 0.0
+#endif
          ENDIF
 
Index: /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/bias.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/bias.F90	(revision 8400)
+++ /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/bias.F90	(revision 8400)
@@ -0,0 +1,866 @@
+MODULE bias 
+   !! Is used by OPA and STEP
+   !!======================================================================
+   !!                 *** Module bias ***
+   !! Code to estimate and apply bias correction.
+   !! The bias is in T/S and Pressure. It follows the generalized
+   !! bias algorithm presented in Balmaseda et al 2007.
+   !!
+   !! It can be read from a file offline, estimated online from relaxation
+   !! terms or from assimilation increments (this latter estimtd in inner loop)
+   !!
+   !! The parameters for the time evolution of the bias can be different
+   !! from the relaxation terms and for the assim increments. Only the 
+   !! parameter for the relaxtion terms are considered here.
+   !! 
+   !! The offline bias term can contain the seasonal cycle.
+   !!
+   !! The time evolution of the bias for relaxtion is estimated as followed
+   !!   bias_rlx(t+1)=t_rlx_mem*bias_rlx(t)+t_rlx_upd*(t/s)trdmp.
+   !!
+   !! The total bias in T/S is partion between the correction to T/S only
+   !!  (fb_t) and the correction applied to the pressure gradient (fb_p).
+   !!  We impose that (fb_p=1.-fb_t). These factors can be different for the
+   !!  different terms(fb_t_rxl,fb_t_asm,fb_t_ofl)
+   !!
+   !!    (t/s)bias =  fb_t_ofl * (t/s)bias_ofl +
+   !!                 fb_t_rlx * (t/s)bias_rlx +
+   !!                 fb_t_asm * (t/s)bias_asm
+   !!    (t/s)bias_p =fb_p_ofl * (t/s)bias_ofl+ 
+   !!                 fb_p_rlx * (t/s)bias_rlx_p +
+   !!                 fb_p_asm * (t/s)bias_asm_p
+   !!    (t/s)bias is applied directely to correct T and S
+   !!    (t/s)bias_p is used to compute rhd_pc and gru/v_pc
+   !!
+   !!  Note: the above is an adhoc /simple way of making the partition
+   !!        between bias in pressure and bias in T/S. It would be better
+   !!        if the partition is done at the time of estimating the time
+   !!        evolution of the bias. That would mean doubling the number of
+   !!        3D arrays.
+   !!
+   !!  New addtion: (not in Balmaseda et al 2007):
+   !!  A more physical alternative to the partition of the bias can be
+   !!  done in terms of the inertial frequency: when the time scales of 
+   !!  adjustment are shorter that >1/f (Equator), the bias correction should
+   !!  be in the the pressure term. Otherwise, it can act directly in T/S.
+   !!  NOTE that the bias correction in the pressure term here (following
+   !!  (Bell et al 2007) is just the "opposite" as the semi-prognostic method 
+   !!  in Greatbatch et al 2004. 
+   !!  The use of this partition is controlled by ln_inertial=.true.
+   !!  
+   !!
+   !!        2009-03              (M.A. Balmaseda ECMWF)
+   !!======================================================================
+
+   !!----------------------------------------------------------------------
+   !!   bias_init  : Read in the bias namelist and the bias arrays
+   !!   tra_bias   : Apply the bias fields on T/S directly
+   !!   dyn_bias   : Compute density correction for dynamic hpg
+   !!   bias_opn   : open bias files for restart capabilities
+   !!   bias_wrt   : write bias fies "     "          "
+   !!----------------------------------------------------------------------
+   !! * Modules used
+   USE par_kind, ONLY: &
+      & wp
+   USE par_oce, ONLY: &
+      & jpi, &
+      & jpj, &
+      & jpk
+   USE dom_oce, ONLY: &
+      & rdt,          &
+      & ln_zps,       &
+      & gphit
+   USE phycst,  ONLY: &
+      & rday,         &
+      & rad
+   USE oce, ONLY: &
+      & tsb, tsn, tsa, &
+      & rhop,          &
+      & gtsu, gtsv
+   USE dynhpg, ONLY:   &
+      & ln_dynhpg_imp      
+   USE tradmp
+   USE dtatsd, ONLY: &
+      & ln_tsd_tradmp
+   USE in_out_manager, ONLY: &
+      & lwp,          &
+      & numnam_ref,   &
+      & numnam_cfg,   &
+      & numond,       &
+      & numout,       &
+      & lrst_oce,     &
+      & nit000
+   USE iom
+   USE eosbn2
+   USE zpshde          ! partial step: hor. derivative (zps_hde routine)
+   USE biaspar
+   USE fldread         ! read input fields
+   USE lbclnk          ! lateral boundary conditions (or mpp link)
+   USE asmpar
+   USE asminc
+   USE lib_mpp, ONLY: &
+      & ctl_stop, &
+      & ctl_nam
+
+   IMPLICIT NONE
+
+   !! * Routine accessibility
+   PRIVATE   
+   PUBLIC bias_init,   &   !: Read in the bias namelist and the bias arrays
+      &   tra_bias,    &   !: Estimate/apply bias on traces
+      &   dyn_bias,    &   !: " density correction for pressure gradient.
+      &   bias_opn,    &
+      &   bias_wrt
+
+   !! * Substitutions needed to have the variable fsdept_n
+#  include "domzgr_substitute.h90"
+   !! * Shared variables
+   !! * Private module variables
+   REAL(wp), PRIVATE :: &
+      & bias_time_unit_asm,   &  !: bias_asm units in s ( per day = 86400 s)   
+      & bias_time_unit_rlx,   &  !: bias_rlx units in s ( 1 second) 
+      & bias_time_unit_ofl,   &  !: bias_ofl units in s ( 1 second) 
+      & t_rlx_mem,            &  !: time param for mem in bias_rlx model
+      & t_rlx_upd,            &  !: time param for update in bias_rlx model
+                                 !: (pct of incr for computation of bias)
+      & t_asm_mem,            &  !: time param for mem in bias_asm model
+      & t_asm_upd,            &  !: time param for update in bias_asm model
+                                 !: (pct of incr for computation of bias)
+      & fb_t_rlx,             &  !: parition of bias in T for rlx bias term
+      & fb_t_asm,             &  !: parition of bias in T for asm bias term
+      & fb_t_ofl,             &  !: parition of bias in T for ofl bias term
+      & fb_p_rlx,             &  !: parition of bias in P for rlx bias term
+      & fb_p_asm,             &  !: parition of bias in P for asm bias term
+      & fb_p_ofl,             &  !: parition of bias in P for ofl bias term
+      & fctamp,               &  !: amplification factor for T if inertial
+      & rn_maxlat_bias,         &  !: Max lat for latitudinal ramp
+      & rn_minlat_bias             !: Min lat for latitudinal ramp
+
+   LOGICAL,  PRIVATE :: lalloc
+   REAL(wp), PRIVATE, DIMENSION(:,:,:), ALLOCATABLE :: &
+      & tbias_asm, &       !: Temperature bias field
+      & sbias_asm, &       !: Salinity bias field
+      & tbias_rlx, &       !: Temperature bias field
+      & sbias_rlx, &       !: Salinity bias field
+      & tbias_asm_out, &   !: Output temperature bias field
+      & sbias_asm_out, &   !: Output salinity bias field
+      & tbias_rlx_out, &   !: Output temperature bias field
+      & sbias_rlx_out, &   !: Output salinity bias field
+      & tbias_p_out,   &   !: Output temperature bias field for P correction
+      & sbias_p_out        !: Output salinity bias field for P correction
+
+   INTEGER, PRIVATE :: nn_lat_ramp     ! choice of latitude dependent ramp
+                                       ! for the pressure correction.
+				       ! 1:inertial ramp, 2:linear ramp, else:no ramp
+   LOGICAL, PRIVATE :: ln_bsyncro      ! syncronous or assincrous bias correction   
+   LOGICAL, PRIVATE :: ln_itdecay      ! evolve bias correction at every time step.   
+
+   REAL(wp), PRIVATE, DIMENSION(:,:), ALLOCATABLE :: fbcoef  
+
+   INTEGER, PRIVATE ::  &
+      & numbias_asm, &    ! netcdf id of bias file from assim
+      & numbias_tot, &    ! netcdf id of bias file with total bias
+      & nn_bias_itwrt     ! time step for outputting bias pressure corr
+         
+   CHARACTER(LEN=128), PRIVATE :: &
+      & cn_bias_asm,   &  ! name of bias file from assim
+      & cn_bias_tot       ! name of bias with total/rlx bias
+
+   ! Structure of input T and S bias offline (file informations, fields read)
+   TYPE(FLD), PRIVATE, ALLOCATABLE, DIMENSION(:) ::   sf_tbias_ofl  
+   TYPE(FLD), PRIVATE, ALLOCATABLE, DIMENSION(:) ::   sf_sbias_ofl
+
+   TYPE(FLD_N), PRIVATE ::&   ! information about the fields to be read
+      &  sn_tbias_ofl, sn_sbias_ofl  
+   
+CONTAINS
+
+   SUBROUTINE bias_init
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE bias_init  ***
+      !!
+      !! ** Purpose : Read in the bias namelist and read in the bias fields.
+      !!
+      !! ** Method  : Read in the bias namelist and read in the bias fields.
+      !!
+      !! ** Action  :
+      !!
+      !! History :
+      !!        !  08-05  (D. Lea)    Initial version
+      !!        !  08-10  (M. Martin) Tidy
+      !!        !  09-03  (M. Balmaseda). Generalize to estimate the bias
+      !!                                  from relax and offline bias term.
+      !!                                  Introduce parameters to control the
+      !!                                  model for the bias 
+      !!                                  (variables and time evolution)
+      !!----------------------------------------------------------------------
+
+      IMPLICIT NONE
+      
+      !! * Local declarations
+
+      CHARACTER(len=100) ::  cn_dir       ! dir for location ofline bias
+      INTEGER            ::  ierror
+      INTEGER            ::  ios          ! Local integer output status for namelist read
+      REAL(wp)           ::  eft_rlx,  &  ! efolding time (bias memory) in days
+         &                   eft_asm,  &  !     "      "
+         &                   log2,     &
+         &                   lenscl_bias     !lengthscale of the pressure bias decay between minlat and maxlat.
+      
+      NAMELIST/nambias/ ln_bias, ln_bias_asm, ln_bias_rlx, ln_bias_ofl,   &
+         & ln_bias_ts_app, ln_bias_pc_app,                                &
+         & fb_t_asm, fb_t_rlx, fb_t_ofl, fb_p_asm, fb_p_rlx, fb_p_ofl,    &
+         & eft_rlx, t_rlx_upd, eft_asm, t_asm_upd, nn_lat_ramp,           &
+         & bias_time_unit_asm, bias_time_unit_rlx, bias_time_unit_ofl,    &
+         & cn_bias_tot, cn_bias_asm, cn_dir, sn_tbias_ofl, sn_sbias_ofl,  &
+         & ln_bsyncro, fctamp, rn_maxlat_bias, rn_minlat_bias,            &
+         & nn_bias_itwrt, ln_itdecay
+ 
+
+      !-----------------------------------------------------------------------
+      ! Read Namelist : bias interface
+      !-----------------------------------------------------------------------
+
+
+      REWIND( numnam_ref )              ! Namelist nambias in reference namelist : Bias pressure correction
+      READ  ( numnam_ref, nambias, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambias in reference namelist', lwp )
+
+
+      ! Set additional default values (note that most values are set in the reference namelist)
+      
+      IF ( ln_asmiau ) nn_bias_itwrt = nitiaufin
+      
+      ! ... default values (NB: frequency positive => hours, negative => months)
+      !            !   file    ! frequency !  variable  ! time intep !  clim   ! 'yearly' or !
+      !            !   name    !  (hours)  !   name     !   (T/F)    !  (T/F)  !  'monthly'  !
+      sn_tbias_ofl = FLD_N( 'tbias_ofl'    ,    -1.    ,  'tbias'     ,  .TRUE.   , .FALSE. ,   'yearly', '', '', ''  )
+      sn_sbias_ofl = FLD_N( 'sbias_ofl'    ,    -1.    ,  'sbias'     ,  .TRUE.   , .FALSE. ,   'yearly', '', '', ''  )
+
+
+      REWIND( numnam_cfg )              ! Namelist nambias in configuration namelist : Bias pressure correction
+      READ  ( numnam_cfg, nambias, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambias in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, nambias )
+
+
+      IF ( ( .NOT. ln_bias_asm ) .AND. ( .NOT. ln_bias_ofl ) .AND. ( .NOT. ln_bias_rlx ) ) THEN
+         ln_bias_ts_app = .FALSE.
+         ln_bias_pc_app = .FALSE.
+         ln_bias        = .FALSE.
+      ENDIF
+      
+      ! set up decay scales
+      log2           = LOG( 2.0_wp )
+      t_rlx_mem      = EXP( - log2 * rdt / ( eft_rlx * rday ) )
+      t_asm_mem      = EXP( - log2 * bias_time_unit_asm/ ( eft_asm * rday ) )
+      
+      ! Control print
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'bias_init : '
+         WRITE(numout,*) '~~~~~~~~~~~ '
+         WRITE(numout,*) '  Namelist nambias : '
+
+         WRITE(numout,*) '  Bias switches/options/variables '
+         WRITE(numout,*) '     bias main switch               ln_bias        = ',ln_bias
+         WRITE(numout,*) '     bias from assim                ln_bias_asm    = ',ln_bias_asm
+         WRITE(numout,*) '     bias from relax                ln_bias_rlx    = ',ln_bias_rlx
+         WRITE(numout,*) '     bias from offln                ln_bias_ofl    = ',ln_bias_ofl
+         WRITE(numout,*) '     bias T and S apply             ln_bias_ts_app = ',ln_bias_ts_app
+         WRITE(numout,*) '     bias pressure correctn apply   ln_bias_pc_app = ',ln_bias_pc_app
+         WRITE(numout,*) '     bias pressure correctn apply   ln_bias_pc_app = ',ln_bias_pc_app
+         WRITE(numout,*) '     lat ramp for bias correction   nn_lat_ramp    = ',nn_lat_ramp
+         WRITE(numout,*) '     time step for writing bias fld nn_bias_itwrt  = ',nn_bias_itwrt
+	 WRITE(numout,*) '     evolve pcbias at each timestep ln_itdecay     = ',ln_itdecay
+         WRITE(numout,*) '  Parameters for parition of bias term '
+         WRITE(numout,*) '                                    fb_t_rlx       = ',fb_t_rlx
+         WRITE(numout,*) '                                    fb_t_asm       = ',fb_t_asm
+         WRITE(numout,*) '                                    fb_t_ofl       = ',fb_t_ofl
+         WRITE(numout,*) '                                    fb_p_rlx       = ',fb_p_rlx
+         WRITE(numout,*) '                                    fb_p_asm       = ',fb_p_asm
+         WRITE(numout,*) '                                    fb_p_ofl       = ',fb_p_ofl
+         WRITE(numout,*) '  Parameters for time evolution of bias '
+         WRITE(numout,*) '  Rlx   efolding time (mem)         eft_rlx,t_rlx_mem = ', eft_rlx, t_rlx_mem, 1. - log2 * rdt / (eft_rlx * rday)
+         WRITE(numout,*) '        uptdate factor              t_rlx_upd         = ',t_rlx_upd
+         WRITE(numout,*) '  Asm   efolding time (mem)         eft_asm,t_asm_mem = ', eft_asm, t_asm_mem, 1. - log2 * rdt / (eft_asm * rday)
+         WRITE(numout,*) '        uptdate factor              t_asm_upd         = ',t_asm_upd
+         WRITE(numout,*) '  Filenames and input structures'
+         WRITE(numout,*) '     bias_tot filename              cn_bias_to     = ',cn_bias_tot
+         WRITE(numout,*) '     bias_asm filename              cn_bias_asm    = ',cn_bias_asm
+         WRITE(numout,*) '     bias_asm time unit (secs)  bias_time_unit_asm = ',bias_time_unit_asm
+         WRITE(numout,*) '     structure Tem bias ofl         sn_tbias_ofl   = ',sn_tbias_ofl
+         WRITE(numout,*) '     structure Sal bias ofl         sn_sbias_ofl   = ',sn_sbias_ofl
+         
+         IF ( ( (.NOT. ln_tsd_tradmp) .OR. (.NOT. ln_tradmp) ) .AND. ln_bias_rlx ) &
+            &   CALL ctl_stop (' lk_dtatem, lk_dtasal and lk_tradmp need to be true with ln_bias_rlx' )
+
+      ENDIF
+      IF( .NOT. ln_bias ) RETURN
+
+      IF( .NOT. lalloc ) THEN
+
+         ALLOCATE( tbias(jpi,jpj,jpk)  , &
+            &      sbias(jpi,jpj,jpk)  , &
+            &      tbias_p(jpi,jpj,jpk), &
+            &      sbias_p(jpi,jpj,jpk), &
+            &      rhd_pc(jpi,jpj,jpk) , &
+            &      gru_pc(jpi,jpj)     , &
+            &      grv_pc(jpi,jpj)       )
+
+         ALLOCATE( fbcoef(jpi,jpj)       )
+         
+         IF( ln_bias_asm ) ALLOCATE(  tbias_asm(jpi,jpj,jpk),     &
+            &                         sbias_asm(jpi,jpj,jpk),     &
+                                      tbias_asm_out(jpi,jpj,jpk), &  
+                                      sbias_asm_out(jpi,jpj,jpk), &  
+                                      tbias_p_out(jpi,jpj,jpk),   & 			      
+                                      sbias_p_out(jpi,jpj,jpk)    )
+
+         IF( ln_bias_rlx ) ALLOCATE(  tbias_rlx(jpi,jpj,jpk),     &
+            &                         sbias_rlx(jpi,jpj,jpk),     &
+                                      tbias_rlx_out(jpi,jpj,jpk), &  
+                                      sbias_rlx_out(jpi,jpj,jpk)  )
+         lalloc = .TRUE.
+
+      ENDIF
+
+      IF( ln_bias_ofl ) THEN      ! set sf_tbias_ofl and sf_sbias_ofl strctrs
+         !
+         ! tbias
+         !
+         ALLOCATE( sf_tbias_ofl(1), STAT=ierror )
+         IF( ierror > 0 ) THEN
+            CALL ctl_stop( 'bias_init: unable to allocate sf_tbias_ofl structure' )   ;    RETURN
+         ENDIF
+         ALLOCATE( sf_tbias_ofl(1)%fnow(jpi,jpj,jpk) )
+         ALLOCATE( sf_tbias_ofl(1)%fdta(jpi,jpj,jpk,2) )
+         
+         ! fill structure with values and control print
+         CALL fld_fill( sf_tbias_ofl, (/ sn_tbias_ofl /), cn_dir, 'bias_init', 'Offline T bias term ', 'nam_tbias_ofl' )
+         !
+         ! salinity bias
+         ! 
+         ALLOCATE( sf_sbias_ofl(1), STAT=ierror )
+         IF( ierror > 0 ) THEN
+            CALL ctl_stop( 'bias_init: unable to allocate sf_sbias_ofl structure' )   ;   RETURN
+         ENDIF
+         ALLOCATE( sf_sbias_ofl(1)%fnow(jpi,jpj,jpk) )
+         ALLOCATE( sf_sbias_ofl(1)%fdta(jpi,jpj,jpk,2) )
+         
+         ! fill structure with values and control print
+         CALL fld_fill( sf_sbias_ofl, (/ sn_sbias_ofl /), cn_dir, 'bias_init', 'Offline S bias term ', 'nam_sbias_ofl' )
+         
+      ENDIF
+
+      ! Read total bias
+      IF ( ln_bias ) THEN
+         tbias(:,:,:)   = 0.0_wp
+         sbias(:,:,:)   = 0.0_wp
+         tbias_p(:,:,:) = 0.0_wp
+         sbias_p(:,:,:) = 0.0_wp
+         gru_pc(:,:)    = 0.0_wp
+         grv_pc(:,:)    = 0.0_wp
+         IF ( ln_bias_rlx ) THEN
+            tbias_rlx(:,:,:) = 0.0_wp
+            sbias_rlx(:,:,:) = 0.0_wp
+         ENDIF
+         IF ( ln_bias_asm ) THEN   !now rlx and asm bias in same file
+           tbias_asm(:,:,:) = 0.0_wp
+           sbias_asm(:,:,:) = 0.0_wp
+         ENDIF
+         numbias_tot    = 0
+         ! Get bias from file and prevent fail if the file does not exist
+         IF(lwp) WRITE(numout,*) 'Opening ',TRIM( cn_bias_tot ) 
+         CALL iom_open( cn_bias_tot, numbias_tot, ldstop=.FALSE. )
+         
+         IF ( numbias_tot > 0 ) THEN      
+            ! Could check validity time of bias fields here...
+            ! Get the T and S bias data
+            IF(lwp) WRITE(numout,*) 'Reading bias fields from tot...'
+
+            !Search for bias from relaxation term if needed. Use same file
+            IF ( ln_bias_rlx ) THEN
+               IF(lwp) WRITE(numout,*) 'Reading bias fields for bias rlx from file ',cn_bias_tot
+               IF( iom_varid( numbias_tot, 'tbias_rlx' ) > 0 ) THEN
+                  ! Get the T and S bias data
+                  CALL iom_get( numbias_tot, jpdom_autoglo, 'tbias_rlx', tbias_rlx )
+                  CALL iom_get( numbias_tot, jpdom_autoglo, 'sbias_rlx', sbias_rlx )
+               ELSE
+                  CALL ctl_stop( 'Bias relaxation variables not found in ',cn_bias_tot )
+               ENDIF
+            ENDIF
+
+
+            !Search for bias from assim term if needed. Use same file
+            IF ( ln_bias_asm .and. .not.ln_bsyncro ) THEN
+               IF(lwp) WRITE(numout,*) 'Reading a-syncro bias fields for bias asm from file ',cn_bias_tot
+               IF( iom_varid( numbias_tot, 'tbias_asm' ) > 0 ) THEN
+                  ! Get the T and S bias data
+                  CALL iom_get( numbias_tot, jpdom_autoglo, 'tbias_asm', tbias_asm )
+                  CALL iom_get( numbias_tot, jpdom_autoglo, 'sbias_asm', sbias_asm )
+               ELSE
+                  CALL ctl_stop( 'Bias assim variables not found in ',cn_bias_tot )
+               ENDIF
+            ENDIF
+
+            ! Close the file
+            CALL iom_close(numbias_tot)
+
+         ELSE 
+            IF(lwp) WRITE(numout,*) 'No bias file found so T and S bias fields are set to zero'
+         ENDIF
+
+      ENDIF
+
+     ! for the time being, the bias_asm is read in the same file as
+     ! bias_rlx
+     ! Implications: bias_asm is estimated/evolved in time in the second outer 
+     !               loop only, when the assimilation increments are ready.
+     !               bias_asm is kept cte during the first outer loop.
+     !              => Assyncronous bias correction.
+     ! Alternative: Syncronous bias correction:
+     !              bias_asm estimated/evolved in the first outer loop
+     !              with the asm incrments of the previous cycle.
+     !              bias_asm kept cte during the second outer loop.
+     !              Implication: bias_asm should be estimated really in the
+     !              inner loop.
+      IF ( ln_bsyncro ) THEN
+      ! Read bias from assimilation from a separate file
+      IF ( ln_bias_asm ) THEN
+         tbias_asm(:,:,:) = 0.0_wp
+         sbias_asm(:,:,:) = 0.0_wp
+         numbias_asm      = 0
+         ! Get bias from file and prevent fail if the file does not exist
+         IF(lwp) WRITE(numout,*) 'Opening file for syncro assim bias ',TRIM( cn_bias_asm ) 
+         CALL iom_open( cn_bias_asm, numbias_asm, ldstop=.FALSE. )
+         
+         IF ( numbias_asm > 0 ) THEN      
+            ! Could check validity time of bias fields here...
+            
+            ! Get the T and S bias data
+            IF(lwp) WRITE(numout,*) 'Reading syncro bias fields from asm from file ',cn_bias_asm
+            CALL iom_get( numbias_asm, jpdom_autoglo, 'tbias_asm', tbias_asm )
+            CALL iom_get( numbias_asm, jpdom_autoglo, 'sbias_asm', sbias_asm )
+           
+!  this is only applicable if tbias_asm were to be calculated in the inner loop
+            tbias_asm(:,:,:) = tbias_asm(:,:,:) * rdt / bias_time_unit_asm
+            sbias_asm(:,:,:) = sbias_asm(:,:,:) * rdt / bias_time_unit_asm
+            
+            ! Close the file
+            CALL iom_close(numbias_asm)
+            
+         ELSE 
+            IF(lwp) WRITE(numout,*) 'No bias file found from asm so T and S bias fields are set to zero'
+         ENDIF
+         
+      ENDIF
+
+      ENDIF
+
+      !latitudinal dependence of partition coeficients. Adhoc
+      IF ( nn_lat_ramp == 1 ) THEN
+         ! Use the inertial ramp.
+         lenscl_bias = ( rn_maxlat_bias - rn_minlat_bias )*2._wp
+         WHERE ( abs( gphit(:,:) ) <= rn_minlat_bias )         
+            fbcoef(:,:) = 0._wp          
+         ELSEWHERE ( abs( gphit(:,:) ) >= rn_maxlat_bias )         
+            fbcoef(:,:) = 1._wp                    
+         ELSEWHERE        
+            fbcoef(:,:) = 1._wp - exp( -( abs( gphit(:,:) ) - rn_minlat_bias ) &
+                           * ( abs( gphit(:,:) ) - rn_minlat_bias ) / lenscl_bias )                         
+         ENDWHERE 
+      ELSEIF ( nn_lat_ramp == 2 ) THEN   
+         ! Use a linear ramp consist with the geostrophic velocity balance ramp in NEMOVAR
+     
+         WHERE ( abs( gphit(:,:) ) <= rn_minlat_bias )
+            fbcoef(:,:) = 0._wp
+         ELSEWHERE ( abs( gphit(:,:) ) >= rn_maxlat_bias ) 
+            fbcoef(:,:) = 1._wp
+         ELSEWHERE
+            fbcoef(:,:) = 1._wp - ((rn_maxlat_bias - abs( gphit(:,:)))/(rn_maxlat_bias - rn_minlat_bias))
+         ENDWHERE
+      ELSE
+         fbcoef(:,:) = 0.0_wp
+         fctamp      = 0.0_wp
+      ENDIF
+
+
+   END SUBROUTINE bias_init
+
+   SUBROUTINE tra_bias ( kt )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE tra_bias  ***
+      !!
+      !! ** Purpose : Update bias field every time step
+      !!
+      !! ** Method  : add contributions to bias from 3 terms
+      !!
+      !! ** Action  : Bias from assimilation (read in bias_init) 
+      !!              Bias from relaxation term is estimated according to
+      !!              the prescribed time evolution of the bias
+      !!              Bias from ofl term is read from external file
+      !!              The difference contributions are added and the partition
+      !!              into direct bias in T/S and pressure perfomed.
+      !!
+      !! History :  09-03  (M. Balmaseda) 
+      !!----------------------------------------------------------------------
+      !! called every timestep after dta_sst if ln_bias true.
+
+      IMPLICIT NONE
+
+      !! * Arguments
+      INTEGER, INTENT(IN) ::   kt             ! ocean time-step index
+      !! * Local variables
+      INTEGER             ::   ji,jj,jk, it   ! local loop index
+      REAL(wp)            ::   tsclf          ! time saling factor 
+      REAL(wp)            ::   fb_t_asm_max, fb_t_rlx_max, fb_t_ofl_max
+      REAL(wp)            ::   ztfrac, ztsday
+      REAL(wp), DIMENSION(jpi,jpj) :: zcof1, zcof2
+
+      IF ( .NOT. ln_bias ) RETURN
+      fb_t_rlx_max   = MIN(fb_t_rlx*fctamp,1.0_wp)
+      fb_t_asm_max   = MIN(fb_t_asm*fctamp,1.0_wp)
+      fb_t_ofl_max   = MIN(fb_t_ofl*fctamp,1.0_wp)
+
+      tbias(:,:,:)   = 0.0_wp
+      sbias(:,:,:)   = 0.0_wp
+      tbias_p(:,:,:) = 0.0_wp
+      sbias_p(:,:,:) = 0.0_wp
+      IF ( ln_bias_asm ) THEN
+      
+         tsclf = 1
+         IF ( .NOT.ln_bsyncro ) tsclf = rdt / bias_time_unit_asm 
+         zcof1(:,:) = tsclf * ( ( 1.0_wp - fbcoef(:,:) ) * fb_t_asm + &
+            &                              fbcoef(:,:)   * fb_t_asm_max )
+         zcof2(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_p_asm
+      
+         IF ( ln_itdecay ) THEN   !decay the pressure correction at each time step
+	 
+	    ztsday  = rday / real(rdt)
+
+            IF( ln_asmiau .and. ln_trainc ) THEN  ! nudge in increments and decay historical contribution
+
+               
+               IF ( kt <= nitiaufin ) THEN  ! During IAU calculate the fraction of increments to apply at each time step
+
+                  ztfrac = real(kt) / real(nitiaufin)  ! nudging factor during the IAU
+
+	       
+                  IF (lwp) THEN
+                     WRITE(numout,*) 'tra_bias : bias weights'
+                     WRITE(numout,*) '~~~~~~~~~~~~'
+                     WRITE(numout,* ) "proportion of  increment applied in pcbias ",ztfrac
+                     WRITE(numout,* ) "proportion of  historical pcbias applied ",t_asm_mem**(real(kt)/ztsday)
+                  ENDIF
+	        
+                  DO jk = 1, jpkm1
+                     tbias(:,:,jk) = tbias(:,:,jk) +                            &
+                     &                ( t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk)  +                    &
+                     &                ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof1(:,:)
+                     sbias(:,:,jk) = sbias(:,:,jk) +                            &
+                     &               ( t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk)  +                     &
+                     &               ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof1(:,:)
+                     tbias_p(:,:,jk) = tbias_p(:,:,jk) +                        &
+                     &               ( t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk)  +                     &
+                     &               ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof2(:,:)
+                     sbias_p(:,:,jk) = sbias_p(:,:,jk) +                        &  
+                     &               ( t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk)  +                     &
+                     &               ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof2(:,:)
+                  ENDDO
+
+                  IF ( .not.ln_bsyncro ) THEN 
+                     IF ( kt == nn_bias_itwrt ) THEN
+                        DO jk = 1, jpkm1
+                           tbias_asm_out(:,:,jk) =  t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk)  +       &
+                           &                     ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk)
+                           sbias_asm_out(:,:,jk) =  t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk)  +       &
+                           &                     ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk)
+                         END DO
+                     ENDIF
+                  ENDIF
+
+                  ! update the historical component with the increments at the end of the IAU
+                  IF ( kt == nitiaufin ) THEN
+                     DO jk = 1, jpkm1
+                        tbias_asm(:,:,jk) =  t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk)  +       &
+                        &                     ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk)
+                        sbias_asm(:,:,jk) =  t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk)  +       &
+                        &                     ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk)
+                     END DO
+                  ENDIF
+               
+               ELSE ! decay pressure correction from combined historical component and increments after IAU
+
+                  ztfrac=t_asm_mem**(real(kt-nitiaufin)/real(nitiaufin))  ! decay from end of IAU
+
+                  DO jk = 1, jpkm1
+                     tbias(:,:,jk) = tbias(:,:,jk) +                            &
+                     &                ( ztfrac * tbias_asm(:,:,jk) ) * zcof1(:,:)
+                     sbias(:,:,jk) = sbias(:,:,jk) +                            &
+                     &               (  ztfrac * sbias_asm(:,:,jk) ) * zcof1(:,:)
+                     tbias_p(:,:,jk) = tbias_p(:,:,jk) +                        &
+                     &               (  ztfrac * tbias_asm(:,:,jk) ) * zcof2(:,:)
+                     sbias_p(:,:,jk) = sbias_p(:,:,jk) +                        &  
+                     &               ( ztfrac * sbias_asm(:,:,jk) ) * zcof2(:,:)
+                  ENDDO
+
+                  IF (lwp) THEN
+                     WRITE(numout,*) 'tra_bias : bias weights'
+                     WRITE(numout,*) '~~~~~~~~~~~~'
+                     WRITE(numout,* ) "proportion of increments + historical pcbias applied ",ztfrac
+                  ENDIF
+
+                  IF ( ln_trainc .and. .not.ln_bsyncro ) THEN 
+                     IF ( kt == nn_bias_itwrt ) THEN
+                        DO jk = 1, jpkm1
+                           tbias_asm_out(:,:,jk) =  ztfrac * tbias_asm(:,:,jk) 
+                           sbias_asm_out(:,:,jk) =  ztfrac * sbias_asm(:,:,jk) 
+                         END DO
+                     ENDIF
+                  ENDIF
+
+               ENDIF
+
+
+            ELSE ! no assimilation increments, simply decay pressure correction (e.g for forecasts)
+
+               DO jk = 1, jpkm1
+                  tbias(:,:,jk) = tbias(:,:,jk) +                                                         &
+                  &               ( t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk) ) * zcof1(:,:)
+                  sbias(:,:,jk) = sbias(:,:,jk) +                                                         &
+                  &               ( t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk) ) * zcof1(:,:)
+                  tbias_p(:,:,jk) = tbias_p(:,:,jk) +                                                     &
+                  &               ( t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk) ) * zcof2(:,:)
+                  sbias_p(:,:,jk) = sbias_p(:,:,jk) +                                                     &  
+                  &               ( t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk) ) * zcof2(:,:)
+               ENDDO
+
+               IF (lwp) THEN
+                  WRITE(numout,*) 'tra_bias : bias weights'
+                  WRITE(numout,*) '~~~~~~~~~~~~'
+                  WRITE(numout,* ) "proportion of historical pcbias applied ",t_asm_mem**(real(kt)/ztsday)
+               ENDIF
+
+            ENDIF
+ 
+         ELSE ! maintain a constant pressure correction  
+
+            DO jk = 1, jpkm1
+               tbias(:,:,jk) = tbias(:,:,jk) + tbias_asm(:,:,jk) * zcof1(:,:)
+               sbias(:,:,jk) = sbias(:,:,jk) + sbias_asm(:,:,jk) * zcof1(:,:)
+               tbias_p(:,:,jk) = tbias_p(:,:,jk) + tbias_asm(:,:,jk) * zcof2(:,:)
+               sbias_p(:,:,jk) = sbias_p(:,:,jk) + sbias_asm(:,:,jk) * zcof2(:,:)
+            END DO     
+
+            IF( ln_asmiau .and. ln_trainc .and. .not.ln_bsyncro ) THEN   
+            ! if last outer loop (ln_asmiau=true and ln_trainc=true). t/sbias_asm
+            ! is updated, only once (end of run) taking into account units.
+               IF ( kt == nn_bias_itwrt ) THEN
+                 IF(lwp) WRITE(numout,*)' estimating asm bias at timestep: ',kt
+                 DO jk = 1, jpkm1
+                   tbias_asm_out(:,:,jk) = t_asm_mem * tbias_asm(:,:,jk)  +             &
+                   &                      t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk)
+                   sbias_asm_out(:,:,jk) = t_asm_mem * sbias_asm(:,:,jk) +             &
+                   &                      t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk)               
+                 END DO
+               ENDIF
+             ENDIF
+  
+         ENDIF
+
+      ENDIF
+
+
+#if   defined key_tradmp 
+      ! Time evolution of bias from relaxation
+      IF ( ln_bias_rlx ) THEN
+         tbias_rlx(:,:,:) = t_rlx_mem * tbias_rlx(:,:,:) + &
+            &               t_rlx_upd * ttrdmp(:,:,:) * rdt
+         sbias_rlx(:,:,:) = t_rlx_mem * sbias_rlx(:,:,:) + &
+            &               t_rlx_upd * strdmp(:,:,:) * rdt
+         zcof1(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_t_rlx +      &
+            &                    fbcoef(:,:)   * fb_t_rlx_max
+         zcof2(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_p_rlx
+         DO jk = 1, jpkm1
+            tbias(:,:,jk)   = tbias(:,:,jk) + tbias_rlx(:,:,jk) * zcof1(:,:)
+            sbias(:,:,jk)   = sbias(:,:,jk) + sbias_rlx(:,:,jk) * zcof1(:,:)
+            tbias_p(:,:,jk) = tbias_p(:,:,jk) + tbias_rlx(:,:,jk) * zcof2(:,:)
+            sbias_p(:,:,jk) = sbias_p(:,:,jk) + sbias_rlx(:,:,jk) * zcof2(:,:)
+         ENDDO
+         IF ( kt == nn_bias_itwrt ) THEN
+            tbias_rlx_out(:,:,:) = tbias_rlx(:,:,:)
+            sbias_rlx_out(:,:,:) = sbias_rlx(:,:,:)
+         ENDIF
+      ENDIF
+#endif
+      ! ofline bias
+      IF ( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*) ' tra_bias: ln_bias_ofl = ',ln_bias_ofl
+         IF(lwp) WRITE(numout,*) ' ~~~~~~~~~'
+      ENDIF
+      IF ( ln_bias_ofl ) THEN
+         IF(lwp) WRITE(numout,*) 'reading offline bias'
+         CALL fld_read( kt, 1, sf_tbias_ofl ) 
+         CALL fld_read( kt, 1, sf_sbias_ofl ) 
+
+         zcof1(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_t_ofl +           &
+            &                    fbcoef(:,:)   * fb_t_ofl_max
+         zcof2(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_p_ofl
+         DO jk = 1, jpkm1
+            tbias(:,:,jk)   = tbias(:,:,jk) + sf_tbias_ofl(1)%fnow(:,:,jk) * zcof1(:,:)
+            sbias(:,:,jk)   = sbias(:,:,jk) + sf_sbias_ofl(1)%fnow(:,:,jk) * zcof1(:,:)
+            tbias_p(:,:,jk) = tbias_p(:,:,jk) + sf_tbias_ofl(1)%fnow(:,:,jk) * zcof2(:,:)
+            sbias_p(:,:,jk) = sbias_p(:,:,jk) + sf_sbias_ofl(1)%fnow(:,:,jk) * zcof2(:,:)
+         ENDDO
+      ENDIF
+
+
+      !apply bias on tracers if needed      
+      IF ( ln_bias_ts_app ) THEN
+         
+         ! Add the bias directely to T/s
+         tsa(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) + tmask(:,:,:) * tbias(:,:,:) / rdt
+         tsa(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) + tmask(:,:,:) * sbias(:,:,:) / rdt
+
+         ! lateral boundary conditions (is this needed?)
+         CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1.0_wp )
+         CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1.0_wp )   
+
+      ENDIF
+
+      IF ( kt == nn_bias_itwrt ) THEN
+         tbias_p_out(:,:,:)   = tbias_p(:,:,:)
+         sbias_p_out(:,:,:)   = sbias_p(:,:,:)
+      ENDIF
+
+   END SUBROUTINE tra_bias
+
+   SUBROUTINE dyn_bias( kt )
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE dyn_bias  ***
+      !!
+      !! ** Purpose :   Computes rhd_pc, gru/v_pc bias corrected
+      !!                for hydrostatic pressure gradient
+      !!                depending on time step (semi-implicit/centered)
+      !!                If partial steps computes bottom pressure gradient.
+      !!                These correction terms will affect only the dynamical
+      !!                component (pressure gradient calculation), but not
+      !!                the isopycnal calculation for the lateral mixing.
+      !! 
+      !! ** Method  :   At this stage of the computation, ta and sa are the 
+      !!                after temperature and salinity. If semi-implicit, these
+      !!                are used to compute rho and bottom pressure gradient.
+      !!                If centered, tb,sb are used instead.
+      !!                If bias key is activated, the temperature,salinity are
+      !!                bias corrected in the calculation of the density fields
+      !!                used in the pressure gradient calculation.
+      !!
+      !!
+      !! ** Action  : - rhd_pc ready. rhop will be overwriten later
+      !!              - if ln_zps, bottom density gradients gru/v_pc ready.
+      !!----------------------------------------------------------------------
+      !!
+      !! * Arguments
+      INTEGER, INTENT(IN) ::   kt    ! ocean time-step index
+      !! * Local variables
+      REAL(wp) :: tsw(jpi,jpj,jpk,jpts)
+      !!
+      !!----------------------------------------------------------------------
+      !
+      ! gtu,gsu,gtv,gsv rhop will be overwritten later in step.
+      ! 
+      IF( ln_dynhpg_imp  ) THEN                             ! semi-implicit hpg
+         tsw(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) - tbias_p(:,:,:)
+         tsw(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) - sbias_p(:,:,:)
+      ELSE
+         tsw(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) - tbias_p(:,:,:)
+         tsw(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) - sbias_p(:,:,:)
+      ENDIF
+
+      CALL eos( tsw, rhd_pc, rhop , fsdept_n(:,:,:))
+      
+      ! is this needed?
+      CALL lbc_lnk( rhd_pc, 'T', 1.0_wp )
+
+      ! Partial steps: now horizontal gradient of t,s,rd 
+      ! at the bottom ocean level
+      IF( ln_zps    )  THEN
+         CALL zps_hde( kt, jpts, tsw, gtsu, gtsv,  &
+            &          rhd_pc, gru_pc , grv_pc  )
+      ENDIF
+
+   END SUBROUTINE dyn_bias
+   
+   SUBROUTINE bias_opn( kt )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE bias_opn  ***
+      !!                     
+      !! ** Purpose :  open bias restart file following the same logic as the
+      !!               standard restarts.
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      INTEGER, INTENT(IN) ::   kt     ! ocean time-step
+      !! * Local variables
+      CHARACTER(LEN=20)   ::   clbkt    ! ocean time-step deine as a character
+      CHARACTER(LEN=50)   ::   clbias_tot   ! total bias restart file name
+      !!----------------------------------------------------------------------
+      !
+      IF( lrst_oce .AND. .NOT.lrst_bias ) THEN       ! create bias file
+         IF( nitend > 999999999 ) THEN   ;   WRITE(clbkt, *       ) nitend
+         ELSE                            ;   WRITE(clbkt, '(i8.8)') nitend
+         ENDIF
+         clbias_tot = TRIM(cexper)//"_"//TRIM(ADJUSTL(clbkt))//"_"//TRIM(cn_bias_tot)
+         IF(lwp) THEN
+            WRITE(numout,*)
+            SELECT CASE ( jprstlib )
+            CASE ( jprstdimg )   ;   WRITE(numout,*) '             open tot bias restart binary file: '//clbias_tot
+            CASE DEFAULT         ;   WRITE(numout,*) '             open tot bias restart NetCDF file: '//clbias_tot
+            END SELECT
+         ENDIF
+         CALL iom_open( clbias_tot, numbias_tot , ldwrt = .TRUE., kiolib = jprstlib )
+         lrst_bias=.TRUE.
+
+      ENDIF
+      !
+   END SUBROUTINE bias_opn
+
+   SUBROUTINE bias_wrt( kt )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE bias_wrt  ***
+      !!                     
+      !! ** Purpose :   Write bias restart fields in the format corresponding to jprstlib
+      !!
+      !! ** Method  :   Write in numbias_tot when kt == nitend in output
+      !!                file, save fields which are necessary for restart.
+      !!
+      !! Changed the timestep for writing to nitend rather than nitrst as this causes a
+      !! problem if the nstock list is used to determine the restart output times and
+      !! means that the bias is not output at all. M. Martin. 08/2011.
+      !! Need to check with Magdalena that this is ok for ECMWF.
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      INTEGER, INTENT(IN) ::   kt   ! ocean time-step
+      !!----------------------------------------------------------------------
+      !                                                                     ! the begining of the run [s]
+
+      IF ( ln_bias_rlx ) THEN
+         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_rlx' , tbias_rlx_out )   
+         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_rlx' , sbias_rlx_out )   
+      ENDIF
+      
+      IF ( ln_bias_asm ) THEN
+         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_asm' , tbias_asm_out )   
+         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_asm' , sbias_asm_out )   
+         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_p'   , tbias_p_out )
+         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_p'   , sbias_p_out )         
+      ENDIF
+      
+      IF( kt == nitend ) THEN
+         CALL iom_close( numbias_tot )     ! close the restart file (only at last time step)
+         lrst_bias = .FALSE.
+      ENDIF
+      !
+   END SUBROUTINE bias_wrt
+
+END MODULE bias
Index: /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/biaspar.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/biaspar.F90	(revision 8400)
+++ /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/biaspar.F90	(revision 8400)
@@ -0,0 +1,40 @@
+MODULE biaspar
+   !! Variables relevant to bias module
+   !!======================================================================
+   !!                 *** Module biaspar ***
+   !!----------------------------------------------------------------------
+   !! * Modules used
+   USE par_kind, ONLY: &
+      & wp
+   USE par_oce, ONLY: &
+      & jpi, &
+      & jpj, &
+      & jpk
+
+   IMPLICIT NONE
+   PUBLIC
+
+   !! * Shared module variables
+   LOGICAL, PUBLIC :: ln_bias        = .FALSE. !: estimate (apply) bias arrays 
+   LOGICAL, PUBLIC :: ln_bias_asm    = .FALSE. !: estimate bias from assim incr
+   LOGICAL, PUBLIC :: ln_bias_rlx    = .FALSE. !: estimate bias from relaxation 
+   LOGICAL, PUBLIC :: ln_bias_ofl    = .FALSE. !: bias estimated offline
+   LOGICAL, PUBLIC :: ln_bias_ts_app = .FALSE. !: estimate (apply) bias arrays 
+   LOGICAL, PUBLIC :: ln_bias_pc_app = .FALSE. !: estimate bias from assim incr
+   LOGICAL, PUBLIC :: lrst_bias      = .FALSE. !: estimate bias from assim incr
+
+   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: &
+      & tbias, &       !: Temperature bias field for T correction
+      & tbias_p, &     !:  "           "    "    "   P correction
+      & sbias, &       !: Salinity bias field    for S correction
+      & sbias_p, &     !:  "           "    "    "   P correction
+      & rhd_pc         !: Press corrtd density from online to use in dyn_hpg
+
+
+   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: &
+      & gru_pc, &      !: Press corrtd bottom pressure gradient (x-dir)
+      & grv_pc         !: Press corrtd bottom pressure gradient (y-dir)
+
+
+
+END MODULE biaspar
Index: /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DIA/diainsitutem.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DIA/diainsitutem.F90	(revision 8400)
+++ /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DIA/diainsitutem.F90	(revision 8400)
@@ -0,0 +1,117 @@
+MODULE insitu_tem
+
+  USE dom_oce  ! ocean space and time domain
+  USE oce,            ONLY: tsn
+  USE par_oce,        ONLY: jpi, jpj, jpk, jpkm1
+  USE lbclnk,         ONLY: lbc_lnk
+  USE lib_mpp         ! MPP library
+
+  IMPLICIT NONE
+  PRIVATE
+
+  PUBLIC theta2t
+
+  !! * Accessibility
+  PUBLIC insitu_tem_alloc          ! routines called by step.F90
+
+  REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: insitu_t
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+
+CONTAINS
+
+   REAL FUNCTION insitu_tem_alloc()
+      !!----------------------------------------------------------------------
+      INTEGER, DIMENSION(2) :: ierr
+      !!----------------------------------------------------------------------
+      !
+      ierr = 0
+      !
+      ALLOCATE( insitu_t(jpi,jpj,jpk), STAT=ierr(1) )
+         !
+      insitu_tem_alloc = MAXVAL(ierr)
+      IF( lk_mpp )   CALL mpp_sum( insitu_tem_alloc )
+      !
+  END FUNCTION insitu_tem_alloc
+
+!-----------------------------------------------------------------------------------
+!
+! Calculate the insitu temperature by integrating the adiabatic lapse rate from zero 
+! to pressure at depth of tracer level. Based on UM subroutine POTTEM and function ATG
+!
+! Initial version: D. Acreman June 2006
+!
+!-----------------------------------------------------------------------------------
+
+  SUBROUTINE theta2t()
+
+    INTEGER, PARAMETER :: num_steps=10                ! number of steps in integration
+    INTEGER            :: step                        ! iteration counter
+    INTEGER            :: ji, jj, jk                  ! loop indices
+    REAL(wp), DIMENSION(jpi,jpj,jpk) :: zP            ! pressure (decibars)
+    REAL(wp), DIMENSION(jpi,jpj,jpk) :: zT            ! temperature at pressure p
+    REAL(wp), DIMENSION(jpi,jpj,jpk) :: zTB           ! temperature at p-dp
+    REAL(wp), DIMENSION(jpi,jpj,jpk) :: zTA           ! temperature at p+dp
+    REAL(wp), DIMENSION(jpi,jpj,jpk) :: zDP           ! pressure step
+    REAL(wp), DIMENSION(jpi,jpj,jpk) :: zSS           ! salinity(PSU) - 35.0
+    REAL(wp), DIMENSION(jpi,jpj,jpk) :: zLAPSE        ! adiabatic lapse rate
+
+!CDIR IEXPAND (ATG)
+
+! Set the pressure interval for the integration. The integration is carried out from 
+! zero (pressure at the surface) to the pressure at the depth of the tracer point. The
+! pressure in decibars is represented by the depth in metres. Pressures are "Oceanographic"
+! pressures equal to absolute pressure minus one atmosphere
+     zDP(:,:,:) = 0.0
+     DO jk = 1, jpkm1
+        DO jj = 1, jpj
+           DO ji = 1, jpi
+              ! These loops expanded for case where fsdept may be 1D
+              zDP(ji,jj,jk) = fsdept(ji,jj,jk) / real(num_steps)
+           END DO
+        END DO
+     END DO
+
+! Salinity at each point
+     zSS(:,:,:)      = tsn(:,:,:,jp_sal) - 35.0 
+
+! Set initial values of temperature and pressure. zT is the temperature at pressure zP, 
+! zTB is the temperature at pressure zP-zdP and zTA is the temperature at pressure zP+zdP
+     zT(:,:,:)  = tsn(:,:,:,jp_tem)
+     zP(:,:,:)  = 0.0         ! Pressure at surface
+     CALL ATG(zP, zT, zSS, zLAPSE)
+     zTB(:,:,:) = zT(:,:,:)  - zLAPSE(:,:,:) * zDP(:,:,:)
+
+     interation: DO step=1, num_steps
+        ! Calculate lapse rate (dT/dP) and hence TA
+        CALL ATG(zP, zT, zSS, zLAPSE)
+        zTA(:,:,:) = zTB(:,:,:) + 2.0 * zLAPSE(:,:,:) * zDP(:,:,:)
+        ! Have calculated TB, T and TA for this pressure, now advance solution by dP
+        zP(:,:,:)  = zP(:,:,:)  + zDP
+        zTB(:,:,:) = zT(:,:,:)
+        zT(:,:,:)  = zTA(:,:,:)
+
+     END DO interation
+
+     insitu_t(:,:,:) = zT(:,:,:) * tmask(:,:,:)
+     CALL lbc_lnk( insitu_t,  'T', 1.0)
+
+   END SUBROUTINE theta2t
+
+   SUBROUTINE ATG(P,T,DS,zLAPSE)
+
+     REAL, INTENT(IN   ) :: P(jpi,jpj,jpk)        ! PRESSURE (DECIBARS)
+     REAL, INTENT(IN   ) :: T(jpi,jpj,jpk)        ! TEMPERATURE (DEG C)
+     REAL, INTENT(IN   ) :: DS(jpi,jpj,jpk)       ! SALINITY (PSU) -35.0
+     REAL, INTENT(  OUT) :: zLAPSE(jpi,jpj,jpk)   ! LAPSE RATE
+
+       zLAPSE = ((( -2.1687E-16*T+1.8676E-14)*T-4.6206E-13)*P &
+           + ((2.7759E-12*T-1.1351E-10)*DS+((-5.4481E-14*T    &
+           + 8.733E-12)*T-6.7795E-10)*T+1.8741E-8))*P         &
+           + (-4.2393E-8*T+1.8932E-6)*DS                      &
+           + ((6.6228E-10*T-6.836E-8)*T+8.5258E-6)*T+3.5803E-5 
+
+   END SUBROUTINE ATG
+
+END MODULE insitu_tem
Index: /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 8399)
+++ /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 8400)
@@ -47,6 +47,6 @@
    USE iom
    USE ioipsl
-   USE dynspg_oce, ONLY: un_adv, vn_adv ! barotropic velocities     
-
+   USE dynspg_oce, ONLY: un_adv, vn_adv ! barotropic velocities    
+   USE insitu_tem, ONLY: insitu_t, theta2t
 #if defined key_lim2
    USE limwri_2 
@@ -164,4 +164,6 @@
       
       CALL iom_put( "toce", tsn(:,:,:,jp_tem) )    ! 3D temperature
+      CALL theta2t ! in-situ temperature conversion
+      CALL iom_put( "tinsitu", insitu_t(:,:,:))    ! in-situ temperature
       CALL iom_put(  "sst", tsn(:,:,1,jp_tem) )    ! surface temperature
       IF ( iom_use("sbt") ) THEN
@@ -202,5 +204,5 @@
          CALL iom_put( "taubot", z2d )           
       ENDIF
-         
+      
       CALL iom_put( "uoce", un(:,:,:)         )    ! 3D i-current
       CALL iom_put(  "ssu", un(:,:,1)         )    ! surface i-current
Index: /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 8399)
+++ /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 8400)
@@ -44,4 +44,5 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
+   USE biaspar         ! bias correction variables
 
    IMPLICIT NONE
@@ -84,4 +85,7 @@
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z_rhd_st  ! tmp density storage for pressure corr
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   z_gru_st  ! tmp ua trends storage for pressure corr
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   z_grv_st  ! tmp va trends storage for pressure corr
       !!----------------------------------------------------------------------
       !
@@ -92,4 +96,20 @@
          ztrdu(:,:,:) = ua(:,:,:)
          ztrdv(:,:,:) = va(:,:,:)
+      ENDIF
+      !
+      IF ( ln_bias .AND. ln_bias_pc_app ) THEN
+
+         !Allocate space for tempory variables
+         ALLOCATE( z_rhd_st(jpi,jpj,jpk), &
+            &      z_gru_st(jpi,jpj),     &
+            &      z_grv_st(jpi,jpj)      )
+
+         z_rhd_st(:,:,:) = rhd(:,:,:)     ! store orig density 
+         rhd(:,:,:)      = rhd_pc(:,:,:)  ! use pressure corrected density
+         z_gru_st(:,:)   = gru(:,:)
+         gru(:,:)        = gru_pc(:,:)
+         z_grv_st(:,:)   = grv(:,:)
+         grv(:,:)        = grv_pc(:,:)
+
       ENDIF
       !
@@ -112,4 +132,15 @@
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' hpg  - Ua: ', mask1=umask,   &
          &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+      !
+      IF ( ln_bias .AND. ln_bias_pc_app )  THEN
+         rhd(:,:,:) = z_rhd_st(:,:,:)     ! restore original density
+         gru(:,:)   = z_gru_st(:,:)
+         grv(:,:)   = z_grv_st(:,:)
+
+         !Deallocate tempory variables
+         DEALLOCATE( z_rhd_st,     &
+            &        z_gru_st,     &
+            &        z_grv_st      )
+      ENDIF
       !
       IF( nn_timing == 1 )  CALL timing_stop('dyn_hpg')
Index: /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 8399)
+++ /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 8400)
@@ -74,5 +74,5 @@
       INTEGER, INTENT(in) ::   kt                      ! time step
       ! 
-      INTEGER             ::   jk                      ! dummy loop indice
+      INTEGER             ::   jk                      ! dummy loop indices
       REAL(wp)            ::   z2dt, z1_rau0           ! local scalars
       !!----------------------------------------------------------------------
@@ -94,4 +94,20 @@
       z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
       IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
+
+
+#if defined key_asminc
+      !                                                ! Include the IAU weighted SSH increment
+      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
+         CALL ssh_asm_inc( kt )
+#if defined key_vvl
+! Don't directly adjust ssh but change hdivn at all levels instead
+! In trasbc also add in the heat and salt content associated with these changes at each level  
+        DO jk = 1, jpkm1                                 
+                 hdivn(:,:,jk) = hdivn(:,:,jk) - ( ssh_iau(:,:) / ( ht_0(:,:) + 1.0 - ssmask(:,:) ) ) * ( e3t_0(:,:,jk) / fse3t_n(:,:,jk) ) * tmask(:,:,jk) 
+        END DO
+      ENDIF
+#endif
+#endif
+
 
       !                                           !------------------------------!
@@ -123,11 +139,4 @@
 #endif
 
-#if defined key_asminc
-      !                                                ! Include the IAU weighted SSH increment
-      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
-         CALL ssh_asm_inc( kt )
-         ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:)
-      ENDIF
-#endif
 
       !                                           !------------------------------!
Index: /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90	(revision 8399)
+++ /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90	(revision 8400)
@@ -110,4 +110,14 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice          ! jpi, jpj
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  a_p, ht_p ! Meltpond fraction and depth
+   
+   !
+   
+   !
+#if defined key_asminc
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ndaice_da          !: NEMO fresh water flux to ocean due to data assim
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   nfresh_da          !: NEMO salt flux to ocean due to data assim
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   nfsalt_da          !: NEMO ice concentration change/second from data assim
+#endif
+      
 #endif
    
@@ -162,4 +172,7 @@
                 ss_iov(jpi,jpj)       , fr_iu(jpi,jpj)        , fr_iv(jpi,jpj)        , &
                 a_i(jpi,jpj,ncat)     , topmelt(jpi,jpj,ncat) , botmelt(jpi,jpj,ncat) , &
+#if defined key_asminc
+                ndaice_da(jpi,jpj)    , nfresh_da(jpi,jpj)    , nfsalt_da(jpi,jpj)    , &
+#endif
                 sstfrz(jpi,jpj)       , STAT= ierr(1) )
    ! Alex West: Allocating tn_ice with 5 categories.  When NEMO is used with CICE, this variable
Index: /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90	(revision 8399)
+++ /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90	(revision 8400)
@@ -56,4 +56,7 @@
                 fresh_ai,fhocn_ai,fswthru_ai,frzmlt,             &
                 flatn_f,fsurfn_f,fcondtopn_f,                    &
+#ifdef key_asminc
+                daice_da,fresh_da,fsalt_da,                    &
+#endif
                 uatm,vatm,wind,fsw,flw,Tair,potT,Qa,rhoa,zlvl,   &
                 swvdr,swvdf,swidr,swidf,Tf,                      &
@@ -301,4 +304,12 @@
  
       CALL wrk_dealloc( jpi,jpj, ztmp1, ztmp2 )
+
+#if defined key_asminc
+	   ! Initialize fresh water and salt fluxes from data assim   
+	   !  and data assimilation index to cice 
+	   nfresh_da(:,:) = 0.0   
+	   nfsalt_da(:,:) = 0.0   
+	   ndaice_da(:,:) = 0.0         
+#endif
       !
       ! In coupled mode get extra fields from CICE for passing back to atmosphere
@@ -454,4 +465,10 @@
       ENDIF
 
+#if defined key_asminc
+!Ice concentration change (from assimilation)
+      ztmp(:,:)=ndaice_da(:,:)*tmask(:,:,1)
+      Call nemo2cice(ztmp,daice_da,'T', 1. )
+#endif 
+
 ! Snowfall
 ! Ensure fsnow is positive (as in CICE routine prepare_forcing)
@@ -716,4 +733,10 @@
          snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) / dt
       ENDIF
+
+#if defined key_asminc
+! Import fresh water and salt flux due to seaice da
+      CALL cice2nemo(fresh_da, nfresh_da,'T',1.0)
+      CALL cice2nemo(fsalt_da, nfsalt_da,'T',1.0)
+#endif
 
 ! Release work space
Index: /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 8399)
+++ /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 8400)
@@ -33,4 +33,7 @@
    USE timing          ! Timing
    USE eosbn2
+#if defined key_asminc   
+   USE asminc          ! Assimilation increment
+#endif
 
    IMPLICIT NONE
@@ -120,4 +123,5 @@
       REAL(wp) ::   zfact, z1_e3t, zdep
       REAL(wp) ::   zalpha, zhk
+      REAL(wp) ::  zt_frz, zpress
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
       !!----------------------------------------------------------------------
@@ -283,4 +287,23 @@
       IF( iom_use('rnf_x_sss') )   CALL iom_put( "rnf_x_sss", rnf*tsn(:,:,1,jp_sal) )   ! runoff term on sss
 
+#if defined key_asminc
+! WARNING: THIS MAY WELL NOT BE REQUIRED - WE DON'T WANT TO CHANGE T&S BUT THIS MAY COMPENSATE ANOTHER TERM...
+! Rate of change in e3t for each level is ssh_iau*e3t_0/ht_0
+! Contribution to tsa should be rate of change in level / per m of ocean? (hence the division by fse3t_n)
+      IF( ln_sshinc ) THEN         ! input of heat and salt due to assimilation
+         DO jj = 2, jpj 
+            DO ji = fs_2, fs_jpim1
+               zdep = ssh_iau(ji,jj) / ( ht_0(ji,jj) + 1.0 - ssmask(ji, jj) )
+               DO jk = 1, jpkm1
+                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   &
+                                        &            + tsn(ji,jj,jk,jp_tem) * zdep * ( e3t_0(ji,jj,jk) / fse3t_n(ji,jj,jk) )
+                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   &
+                                        &            + tsn(ji,jj,jk,jp_sal) * zdep * ( e3t_0(ji,jj,jk) / fse3t_n(ji,jj,jk) )
+               END DO
+            END DO  
+         END DO  
+      ENDIF
+#endif
+ 
       IF( l_trdtra )   THEN                      ! send trends for further diagnostics
          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
Index: /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 8399)
+++ /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 8400)
@@ -227,4 +227,7 @@
 #endif
       !
+      ! Met Office addition: if failed, return non-zero exit code
+      IF( nstop /= 0 )  CALL exit( 9 ) 
+      !
    END SUBROUTINE nemo_gcm
 
@@ -480,4 +483,5 @@
                             CALL dia_hsb_init   ! heat content, salt content and volume budgets
                             CALL     trd_init   ! Mixed-layer/Vorticity/Integral constraints trends
+                            CALL     bias_init  ! Pressure correction bias
       IF( lk_diaobs     ) THEN                  ! Observation & model comparison
                             CALL dia_obs_init            ! Initialize observational data
@@ -646,4 +650,5 @@
       !!----------------------------------------------------------------------
       USE diawri    , ONLY: dia_wri_alloc
+      USE insitu_tem, ONLY: insitu_tem_alloc
       USE dom_oce   , ONLY: dom_oce_alloc
       USE ldfdyn_oce, ONLY: ldfdyn_oce_alloc
@@ -662,4 +667,5 @@
       ierr =        oce_alloc       ()          ! ocean
       ierr = ierr + dia_wri_alloc   ()
+      ierr = ierr + insitu_tem_alloc()
       ierr = ierr + dom_oce_alloc   ()          ! ocean domain
       ierr = ierr + ldfdyn_oce_alloc()          ! ocean lateral  physics : dynamics
Index: /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 8399)
+++ /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 8400)
@@ -103,4 +103,6 @@
       IF( ln_crs     )       CALL iom_setkt( kstp - nit000 + 1, TRIM(cxios_context)//"_crs" )   ! tell iom we are at time step kstp
 
+      IF( ln_bias )          CALL bias_opn( kstp )
+
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Update data, open boundaries, surface boundary condition (including sea-ice)
@@ -267,4 +269,5 @@
       IF( lk_trabbl      )   CALL tra_bbl    ( kstp )       ! advective (and/or diffusive) bottom boundary layer scheme
       IF( ln_tradmp      )   CALL tra_dmp    ( kstp )       ! internal damping trends
+      IF( ln_bias        )   CALL tra_bias   ( kstp )
       IF( lk_bdy         )   CALL bdy_tra_dmp( kstp )       ! bdy damping trends
                              CALL tra_adv    ( kstp )       ! horizontal & vertical advection
@@ -290,4 +293,5 @@
                &                                           rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
                &                                    gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level
+            IF( ln_bias )    CALL dyn_bias( kstp )
       ELSE                                                  ! centered hpg  (eos then time stepping)
          IF ( .NOT. lk_dynspg_ts ) THEN                     ! eos already called in time-split case
@@ -303,4 +307,5 @@
          IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection
                              CALL tra_nxt( kstp )                ! tracer fields at next time step
+         IF( ln_bias )       CALL dyn_bias( kstp )
       ENDIF
 
@@ -377,4 +382,7 @@
       ENDIF
 
+
+      IF( lrst_bias )          CALL bias_wrt     ( kstp )
+
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Coupled mode
Index: /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 8399)
+++ /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 8400)
@@ -100,5 +100,7 @@
 
    USE crsfld           ! Standard output on coarse grid   (crs_fld routine)
-
+   USE biaspar          ! bias param
+   USE bias             ! bias routines                    (tra_bias routine
+                        !                                   dyn_bias routine)
    USE asminc           ! assimilation increments      (tra_asm_inc routine)
    !                                                   (dyn_asm_inc routine)
