Index: /branches/UKMO/dev_r5518_GO6_package_FOAMv14_updated_pcbias/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_FOAMv14_updated_pcbias/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 11554)
+++ /branches/UKMO/dev_r5518_GO6_package_FOAMv14_updated_pcbias/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 11555)
@@ -85,11 +85,10 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
-      INTEGER                                   ::   ji, jj, jk   ! dummy loop indices
-      INTEGER                                   ::   iku, ikv          ! k indices for bottom level at u and v points
+      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_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
@@ -108,16 +107,21 @@
       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_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(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(:,:,:)
+         z_va_bpc(:,:,:)     = va(:,:,:)
+         
       END IF 
 
@@ -134,18 +138,29 @@
       IF ( ln_bias .AND. ln_bias_pc_app ) THEN
 
-         z_rhd_st(:,:,:) = rhd(:,:,:)     ! store orig density 
-         rhd(:,:,:)      = rhd_pc(:,:,:)  ! use pressure corrected density
+         ! 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(:,:)
-         gru(:,:)        = gru_pc(:,:)
          z_grv_st(:,:)   = grv(:,:)
-         grv(:,:)        = grv_pc(:,:)
-
+
+         ! 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 
+         z_va(:,:,:)     = va(:,:,:)
+
+         ! reset the acceleration trends to their original values
          ua(:,:,:)       = z_ua_bpc(:,:,:)    
-	 va(:,:,:)       = z_va_bpc(:,:,:)  
+         va(:,:,:)       = z_va_bpc(:,:,:)  
 
          ! re-calculate the horizontal pressure gradients with the bpc fields 
@@ -174,16 +189,13 @@
 
          ! 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 
+         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(:,:,:)
-
-         IF(lwp) THEN 
-           WRITE(numout,*) " ! restore original density"
-         ENDIF
 
          ! restore original density, gru and grv fields 
@@ -192,8 +204,8 @@
          grv(:,:)   = z_grv_st(:,:)
 
-         !Deallocate tempory variables
+         ! 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        &
+            &        z_ua_bpc_bot, z_va_bpc_bot        &
             &                                          )
 	    
