Index: NEMO/branches/UKMO/NEMO_4.0.4_FOAM_pcbias/cfgs/SHARED/namelist_ref
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.4_FOAM_pcbias/cfgs/SHARED/namelist_ref	(revision 15601)
+++ NEMO/branches/UKMO/NEMO_4.0.4_FOAM_pcbias/cfgs/SHARED/namelist_ref	(revision 15602)
@@ -1379,2 +1379,36 @@
    cn_storst_out = "restart_sto" !  suffix of stochastic parameter restart file (output)
 /
+!-----------------------------------------------------------------------
+&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.
+   ln_incpc       = .FALSE.
+/
Index: NEMO/branches/UKMO/NEMO_4.0.4_FOAM_pcbias/src/OCE/ASM/bias.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.4_FOAM_pcbias/src/OCE/ASM/bias.F90	(revision 15601)
+++ NEMO/branches/UKMO/NEMO_4.0.4_FOAM_pcbias/src/OCE/ASM/bias.F90	(revision 15602)
@@ -31,5 +31,5 @@
    !!                 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
+   !!    (t/s)bias_p is used to compute rhd_pc
    !!
    !!  Note: the above is an adhoc /simple way of making the partition
@@ -76,6 +76,5 @@
    USE oce, ONLY: &
       & tsb, tsn, tsa, &
-      & rhop,  &
-      & gtsu, gtsv   
+      & rhop 
    USE tradmp
    USE dtatsd, ONLY: &
@@ -326,7 +325,5 @@
             &      tbias_i(jpi,jpj,jpk), &
             &      sbias_i(jpi,jpj,jpk), &
-            &      rhd_pc(jpi,jpj,jpk) , &
-            &      gru_pc(jpi,jpj)     , &
-            &      grv_pc(jpi,jpj)       )
+            &      rhd_pc(jpi,jpj,jpk)   )
 
          ALLOCATE( fbcoef(jpi,jpj), fbcoef_stscale(jpi,jpj) )
@@ -391,6 +388,5 @@
          tbias_i(:,:,:) = 0.0_wp
          sbias_i(:,:,:) = 0.0_wp
-         gru_pc(:,:)    = 0.0_wp
-         grv_pc(:,:)    = 0.0_wp
+	 rhd_pc(:,:,:)  = 0.0_wp
          
          IF ( ln_bias_rlx ) THEN
@@ -896,5 +892,5 @@
       !!                   ***  ROUTINE dyn_bias  ***
       !!
-      !! ** Purpose :   Computes rhd_pc, gru/v_pc bias corrected
+      !! ** Purpose :   Computes rhd_pc bias corrected
       !!                for hydrostatic pressure gradient
       !!                depending on time step (semi-implicit/centered)
@@ -914,5 +910,4 @@
       !!
       !! ** Action  : - rhd_pc ready. rhop will be overwriten later
-      !!              - if ln_zps, bottom density gradients gru/v_pc ready.
       !!----------------------------------------------------------------------
       !!
@@ -926,13 +921,13 @@
       ! gtu,gsu,gtv,gsv rhop will be overwritten later in step.
       ! 
-      tsw(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) - tbias_p(:,:,:)
-      tsw(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) - sbias_p(:,:,:)
+      tsw(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) - tbias_p(:,:,:)
+      tsw(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) - sbias_p(:,:,:)
       IF ( ln_incpc ) THEN
-         tsw(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) - tbias_p(:,:,:) - tbias_i(:,:,:)
-         tsw(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) - sbias_p(:,:,:) - sbias_i(:,:,:)
+         tsw(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) - tbias_p(:,:,:) - tbias_i(:,:,:)
+         tsw(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) - sbias_p(:,:,:) - sbias_i(:,:,:)
          IF ( ln_incpc_only ) THEN
             IF(lwp) WRITE(numout,*) 'ln_incpc_only =', ln_incpc_only, 'tsw updated with IPC only'
-            tsw(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) - tbias_i(:,:,:)
-            tsw(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) - sbias_i(:,:,:)
+            tsw(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) - tbias_i(:,:,:)
+            tsw(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) - sbias_i(:,:,:)
          ENDIF
       ENDIF
@@ -941,15 +936,5 @@
       IF(lwp) WRITE(numout,*) 'dyn_bias( kt ) calculating rhd_pc, kt =', kt
       CALL eos( tsw, rhd_pc, rhop, gdept_n(:,:,:) )
-      
-      CALL lbc_lnk( 'bias', rhd_pc, 'T', 1. )
-
-
-      ! 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
Index: NEMO/branches/UKMO/NEMO_4.0.4_FOAM_pcbias/src/OCE/ASM/biaspar.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.4_FOAM_pcbias/src/OCE/ASM/biaspar.F90	(revision 15601)
+++ NEMO/branches/UKMO/NEMO_4.0.4_FOAM_pcbias/src/OCE/ASM/biaspar.F90	(revision 15602)
@@ -34,9 +34,3 @@
 
 
-   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: NEMO/branches/UKMO/NEMO_4.0.4_FOAM_pcbias/src/OCE/DYN/dynhpg.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.4_FOAM_pcbias/src/OCE/DYN/dynhpg.F90	(revision 15601)
+++ NEMO/branches/UKMO/NEMO_4.0.4_FOAM_pcbias/src/OCE/DYN/dynhpg.F90	(revision 15602)
@@ -46,4 +46,5 @@
    USE timing          ! Timing
    USE iom
+   USE biaspar         ! bias correction variables
 
    IMPLICIT NONE
@@ -92,5 +93,12 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv
+      INTEGER                                   ::   ji, jj, jk                 ! dummy loop indices
+      INTEGER                                   ::   iku, ikv                   ! k indices for bottom level at u and v points
+      INTEGER                                   ::   numf
+      REAL(wp), POINTER, DIMENSION(:,:,:)       ::   ztrdu, ztrdv
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z_rhd_st, z_rhd_diff       ! tmp density storage for pressure corr
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z_ua, z_va                 ! tmp store for ua and va including hpg but not pressure correction 
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z_ua_bpc, z_va_bpc         ! ua calculated with bias pressure correction  
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   z_ua_bpc_bot, z_va_bpc_bot ! bias pc fields calculated at the ocean bottom
       !!----------------------------------------------------------------------
       !
@@ -102,4 +110,24 @@
          ztrdv(:,:,:) = va(:,:,:)
       ENDIF
+      !
+      IF ( ln_bias .AND. ln_bias_pc_app ) THEN
+
+         ! allocate space for tempory variables for the bias pressure correction (bpc)
+         ALLOCATE( z_rhd_st(jpi,jpj,jpk), &
+            &      z_ua(jpi,jpj,jpk),     &
+            &      z_va(jpi,jpj,jpk),     &  
+            &      z_ua_bpc(jpi,jpj,jpk), &
+            &      z_va_bpc(jpi,jpj,jpk), & 
+            &      z_ua_bpc_bot(jpi,jpj), &
+            &      z_va_bpc_bot(jpi,jpj), & 
+	    &      z_rhd_diff(jpi,jpj,jpk)&
+            &    )
+
+         ! save the original acceleration trends 
+         ! (z_ua_bpc, z_va_bpc are used as temporary storage)
+         z_ua_bpc(:,:,:)     = ua(:,:,:)
+         z_va_bpc(:,:,:)     = va(:,:,:)
+         
+      END IF
       !
       SELECT CASE ( nhpg )      ! Hydrostatic pressure gradient computation
@@ -112,4 +140,75 @@
       END SELECT
       !
+      IF ( ln_bias .AND. ln_bias_pc_app ) THEN
+
+         ! The aim here is to calculate the contribution of the bpc to the acceleration terms.
+         ! This is done so that the effect of the bpc on the hpg at the bottom can be removed.
+         ! In order to do that:
+         !    1. The hpg calculation is done again, but with the contributions of the bpc included.
+         !    2. The difference between the acceleration terms (w and w/o bpc) is then calculated.
+         !    3. The effect of the bpc on the bottom hpg is then removed.
+         !    4. The total change to the acceleration terms is then calculated.
+         
+         ! The original density field (without the bpc) is stored.          
+         z_rhd_st(:,:,:) = rhd(:,:,:)
+
+         ! Set the density used in the hpc calculations to the value including the effect of the bpc.
+         rhd(:,:,:)      = rhd_pc(:,:,:) 
+            
+         ! save the acceleration trends including hpg field but calculated without the bpc fields
+         z_ua(:,:,:)     = ua(:,:,:)
+         z_va(:,:,:)     = va(:,:,:)
+
+         ! reset the acceleration trends to their original values
+         ua(:,:,:)       = z_ua_bpc(:,:,:)    
+         va(:,:,:)       = z_va_bpc(:,:,:)  
+
+         ! re-calculate the horizontal pressure gradients with the bpc fields 
+         SELECT CASE ( nhpg )      ! Hydrostatic pressure gradient computation
+         CASE (  np_zco )   ;   CALL hpg_zco    ( kt )      ! z-coordinate
+         CASE (  np_zps )   ;   CALL hpg_zps    ( kt )      ! z-coordinate plus partial steps (interpolation)
+         CASE (  np_sco )   ;   CALL hpg_sco    ( kt )      ! s-coordinate (standard jacobian formulation)
+         CASE (  np_djc )   ;   CALL hpg_djc    ( kt )      ! s-coordinate (Density Jacobian with Cubic polynomial)
+         CASE (  np_prj )   ;   CALL hpg_prj    ( kt )      ! s-coordinate (Pressure Jacobian scheme)
+         CASE (  np_isf )   ;   CALL hpg_isf    ( kt )      ! s-coordinate similar to sco modify for ice shelf
+         END SELECT
+
+         ! calculate the bpc contribution to ua and va
+         z_ua_bpc(:,:,:) = ua(:,:,:) - z_ua(:,:,:)
+         z_va_bpc(:,:,:) = va(:,:,:) - z_va(:,:,:)
+
+         ! calculate the bpc contribution to ua and va at the bottom 
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               iku = mbku(ji,jj)
+               ikv = mbkv(ji,jj) 
+               z_ua_bpc_bot(ji,jj) = z_ua_bpc(ji,jj,iku)
+               z_va_bpc_bot(ji,jj) = z_va_bpc(ji,jj,ikv)
+            END DO  ! ji
+         END DO ! jj
+
+         ! subtract off the bottom values of bpc contribution to ua and va	
+         DO jk = 1, jpk - 1
+            z_ua_bpc(:,:,jk) = z_ua_bpc(:,:,jk) - z_ua_bpc_bot(:,:)
+            z_va_bpc(:,:,jk) = z_va_bpc(:,:,jk) - z_va_bpc_bot(:,:)
+         END DO 
+
+         ! calculate ua using the original hpg (z_ua) and the bias hpg 
+         ! with the bottom pressure gradients subtracted off 
+         ua(:,:,:) = z_ua(:,:,:) + z_ua_bpc(:,:,:)
+         va(:,:,:) = z_va(:,:,:) + z_va_bpc(:,:,:)
+
+         ! restore original density field
+         rhd(:,:,:) = z_rhd_st(:,:,:)     
+
+         ! deallocate tempory variables
+         DEALLOCATE( z_rhd_st, z_ua, z_va, z_ua_bpc,   &
+            &        z_va_bpc, z_ua_bpc_bot,           &
+            &        z_va_bpc_bot, z_rhd_diff          &
+            &                                          )
+	    
+      ENDIF  ! ln_bias .AND. ln_bias_pc_app
+
+
       IF( l_trddyn ) THEN      ! save the hydrostatic pressure gradient trends for momentum trend diagnostics
          ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
Index: NEMO/branches/UKMO/NEMO_4.0.4_FOAM_pcbias/src/OCE/step.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.4_FOAM_pcbias/src/OCE/step.F90	(revision 15601)
+++ NEMO/branches/UKMO/NEMO_4.0.4_FOAM_pcbias/src/OCE/step.F90	(revision 15602)
@@ -179,4 +179,6 @@
                          CALL dyn_ldf       ( kstp )  ! lateral mixing
       IF( ln_zdfosm  )   CALL dyn_osm       ( kstp )  ! OSMOSIS non-local velocity fluxes
+      IF( ln_bias    )   CALL tra_bias      ( kstp )  ! tracer bias
+      IF( ln_bias    )   CALL dyn_bias      ( kstp )
                          CALL dyn_hpg       ( kstp )  ! horizontal gradient of Hydrostatic pressure
                          CALL dyn_spg       ( kstp )  ! surface pressure gradient
@@ -232,5 +234,4 @@
       IF( ln_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 )  ! tracer bias
       IF( ln_bdy     )   CALL bdy_tra_dmp   ( kstp )  ! bdy damping trends
 #if defined key_agrif
@@ -265,5 +266,4 @@
 !!jc2: dynnxt must be the latest call. e3t_b are indeed updated in that routine
                          CALL tra_nxt       ( kstp )  ! finalize (bcs) tracer fields at next time step and swap
-      IF( ln_bias    )   CALL dyn_bias      ( kstp )
                          CALL dyn_nxt       ( kstp )  ! finalize (bcs) velocities at next time step and swap (always called after tra_nxt)
                          CALL ssh_swp       ( kstp )  ! swap of sea surface height
