Index: /branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 11555)
+++ /branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 11556)
@@ -17,4 +17,5 @@
    !!                 !           (A. Coward) suppression of hel, wdj and rot options
    !!            3.6  !  2014-11  (P. Mathiot) hpg_isf: original code for ice shelf cavity
+   !!            3.6  !  2019-08  (M. Bell) Revisions to bias pressure correction 
    !!----------------------------------------------------------------------
 
@@ -84,8 +85,13 @@
       !!----------------------------------------------------------------------
       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
+      INTEGER                                   ::   ji, jj, jk                 ! dummy loop indices
+      INTEGER                                   ::   iku, ikv                   ! k indices for bottom level at u and v points
+      REAL(wp), POINTER, DIMENSION(:,:,:)       ::   ztrdu, ztrdv
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z_rhd_st                   ! 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_gru_st, z_grv_st         ! tmp ua and va trends storage for pressure corr
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   z_ua_bpc_bot, z_va_bpc_bot ! bias pc fields calculated at the ocean bottom
+      
       !!----------------------------------------------------------------------
       !
@@ -98,19 +104,25 @@
       ENDIF
       !
+
       IF ( ln_bias .AND. ln_bias_pc_app ) THEN
 
-         !Allocate space for tempory variables
+         ! allocate space for tempory variables for the bias pressure correction (bpc)
          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
+            &      z_grv_st(jpi,jpj),     &
+            &      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)  & 
+            &    )
+
+         ! 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
@@ -123,4 +135,82 @@
       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 fields etc (without the bpc) are stored.          
+         z_rhd_st(:,:,:) = rhd(:,:,:)
+         z_gru_st(:,:)   = gru(:,:)
+         z_grv_st(:,:)   = grv(:,:)
+
+         ! Set the density etc used in the hpc calculations to the value including the effect of the bpc.
+         rhd(:,:,:)      = rhd_pc(:,:,:)
+         gru(:,:)        = gru_pc(:,:)   
+         grv(:,:)        = grv_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 (  0 )   ;   CALL hpg_zco    ( kt )      ! z-coordinate
+         CASE (  1 )   ;   CALL hpg_zps    ( kt )      ! z-coordinate plus partial steps (interpolation)
+         CASE (  2 )   ;   CALL hpg_sco    ( kt )      ! s-coordinate (standard jacobian formulation)
+         CASE (  3 )   ;   CALL hpg_djc    ( kt )      ! s-coordinate (Density Jacobian with Cubic polynomial)
+         CASE (  4 )   ;   CALL hpg_prj    ( kt )      ! s-coordinate (Pressure Jacobian scheme)
+         CASE (  5 )   ;   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, gru and grv fields 
+         rhd(:,:,:) = z_rhd_st(:,:,:)     
+         gru(:,:)   = z_gru_st(:,:)
+         grv(:,:)   = z_grv_st(:,:)
+
+         ! deallocate tempory variables
+         DEALLOCATE( z_rhd_st,   z_gru_st, z_grv_st,   &
+            &        z_ua, z_va, z_ua_bpc, z_va_bpc,   &
+            &        z_ua_bpc_bot, z_va_bpc_bot        &
+            &                                          )
+	    
+      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(:,:,:)
@@ -132,18 +222,4 @@
       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
-         IF(lwp) THEN 
-         WRITE(numout,*) " ! restore original density"
-         ENDIF
-         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')
