Index: NEMO/branches/UKMO/NEMO_4.0-TRUNK_r14960_HPG/src/OCE/DYN/dynhpg.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0-TRUNK_r14960_HPG/src/OCE/DYN/dynhpg.F90	(revision 15727)
+++ NEMO/branches/UKMO/NEMO_4.0-TRUNK_r14960_HPG/src/OCE/DYN/dynhpg.F90	(revision 15728)
@@ -2,5 +2,6 @@
 
    !!  v2:  takes into account the stretching of the vertical grid in hpg_ffr
-
+   !!  v3:  takes into account the stretching of the vertical grid in the calculation of the reference profile
+   
    !!======================================================================
    !!                       ***  MODULE  dynhpg  ***
@@ -94,7 +95,9 @@
    REAL(wp)    ::   bco_bc_z_srf, aco_bc_z_bot, bco_bc_z_bot                    !    "
 
-   LOGICAL     ::   ln_hpg_djr_ref_ccs  ! T => constrained cubic, F => simple cubic used for interpolation of reference to target profiles
-   LOGICAL     ::   ln_hpg_ffr_ref      ! T => reference profile is subtracted; F => no reference profile subtracted
-   LOGICAL     ::   ln_hpg_ffr_ref_ccs  ! T => use constrained cubic spline to vertically interpolate reference to each profile   
+   LOGICAL     ::   ln_hpg_ref      ! T => reference profile is subtracted; F => no reference profile subtracted
+   LOGICAL     ::   ln_hpg_ref_str  ! T => reference profile calculated taking into account vertical variation in grid spacing
+   LOGICAL     ::   ln_hpg_ref_ccs  ! T => constrained cubic, F => simple cubic used for interpolation of reference to target profiles   
+   LOGICAL     ::   ln_hpg_ref_off  ! only applies if ln_hpg_ref_ccs = T; T => use off-centred simple cubic at upper & lower boundaries 
+
    LOGICAL     ::   ln_hpg_ffr_hor_ccs  ! T => use constrained cubic spline to interpolate z_rhd_pmr along upper faces of cells
    LOGICAL     ::   ln_hpg_ffr_hor_cub  ! T => use simple cubic polynomial to interpolate z_rhd_pmr along upper faces of cells
@@ -195,7 +198,8 @@
          &                 ln_hpg_bcvN_rhd_hor,    ln_hpg_bcvN_rhd_srf,      & 
          &                 ln_hpg_bcvN_rhd_bot,    ln_hpg_bcvN_z_hor,        &                 
-         &                 ln_hpg_djr_ref_ccs,     ln_hpg_ffr_vrt_quad,      &   
-         &                 ln_hpg_ffr_ref,         ln_hpg_ffr_ref_ccs,       & 
+         &                 ln_hpg_ref,             ln_hpg_ref_str,           &
+         &                 ln_hpg_ref_ccs,         ln_hpg_ref_off,           &
          &                 ln_hpg_ffr_hor_ccs,     ln_hpg_ffr_hor_cub,       &
+         &                 ln_hpg_ffr_vrt_quad,                              &   
          &                 ln_dbg_hpg, ln_dbg_ij,  ln_dbg_ik,  ln_dbg_jk,    &
          &                 ki_dbg_min, ki_dbg_max, ki_dbg_cen,               &  
@@ -322,21 +326,22 @@
       END IF
 
-      IF ( lwp .AND. ln_hpg_djr ) THEN
-         IF ( ln_hpg_djr_ref_ccs ) THEN 
-	    WRITE(numout,*) '           interpolation of reference profile by constrained cubic spline'                      
-         ELSE 
+      IF ( lwp .AND. ln_hpg_ref ) THEN
+         IF ( ln_hpg_ref_ccs .AND. ln_hpg_ref_off .AND. ln_hpg_ref_str ) THEN 
+	    WRITE(numout,*) '           interpolation of reference profile by constrained cubic spline with off-centring at boundaries'                  
+         ELSE IF ( ln_hpg_ref_ccs ) THEN
+	    WRITE(numout,*) '           interpolation of reference profile by constrained cubic spline using boundary conditions'
+	 ELSE 
 	    WRITE(numout,*) '           interpolation of reference profile by simple cubic spline with off-centring at boundaries'  
          END IF 
+         IF ( ln_hpg_ref_str ) THEN 
+	    WRITE(numout,*) '           interpolation of reference profile takes account of variation in vertical grid spacing'                  
+	 ELSE 
+	    WRITE(numout,*) '           interpolation of reference profile does NOT take account of variation in vertical grid spacing' 
+         END IF 
       END IF       
 
-      IF ( lwp .AND. ln_hpg_ffr ) THEN
-         IF ( ln_hpg_ffr_ref_ccs ) THEN 
-            IF ( ln_hpg_ffr_ref_ccs ) THEN 
-	       WRITE(numout,*) '           interpolation of reference profile in vertical by constrained cubic spline '
-            ELSE 
-	       WRITE(numout,*) '           interpolation of reference profile in vertical by pure cubic polynomial fit '                   
-            END IF       
-         ELSE 
-	    WRITE(numout,*) '           no reference profile subtracted '                   
+      IF ( lwp .AND. ln_hpg_ffr ) THEN 
+         IF ( .NOT. ln_hpg_ref  ) THEN
+            WRITE(numout,*) '           no reference profile subtracted '                   
          END IF       
          IF ( ln_hpg_ffr_hor_ccs ) THEN 
@@ -1751,12 +1756,12 @@
 	 
       !--------------------------------------------------------------------------------------------------------
-      !  2.2  IF  ln_hpg_djr_ref_ccs compute zdrhd_k_ref then set bcs at top & bottom  
+      !  2.2  IF  ln_hpg_ref_ccs compute zdrhd_k_ref then set bcs at top & bottom  
       !                              (bcs not needed for simple cubic off-centred at boundaries)
       !--------------------------------------------------------------------------------------------------------
 
-         IF ( ln_hpg_djr_ref_ccs ) THEN
+         IF ( ln_hpg_ref_ccs ) THEN
             CALL calc_drhd_k(zrhd_ref, jk_bot_ref, zdrhd_k_ref) 
             IF ( ln_dbg_hpg ) CALL dbg_3dr( '2.3 zdrhd_k_ref', zdrhd_k_ref )	    
-         END IF  ! ln_hpg_djr_ref_ccs 
+         END IF  ! ln_hpg_ref_ccs 
 
       !----------------------------------------------------------------------------------------
@@ -1773,10 +1778,25 @@
             END IF 
 
-            IF ( ln_hpg_djr_ref_ccs ) THEN 
-               CALL ref_to_tgt_ccs ( jio, jjo, gde3w, rhd, zz_ref, zrhd_ref, zdrhd_k_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) )
-            ELSE 
-               CALL ref_to_tgt_cub ( jio, jjo, gde3w, rhd, zz_ref, zrhd_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) )  
+
+            IF ( ln_hpg_ref ) THEN 
+               IF ( ln_hpg_ref_str ) THEN 
+                  IF ( ln_hpg_ref_ccs ) THEN 
+                     IF ( ln_hpg_ref_off ) THEN 
+                        CALL ref_to_tgt_ccs_str_off ( jio, jjo, gde3w, rhd, zz_ref, zrhd_ref, zdrhd_k_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) )
+                     ELSE
+                        CALL ref_to_tgt_ccs_str ( jio, jjo, gde3w, rhd, zz_ref, zrhd_ref, zdrhd_k_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) )
+                     END IF 
+                  ELSE  
+                     CALL ref_to_tgt_cub_str ( jio, jjo, gde3w, rhd, zz_ref, zrhd_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) )  
+                  END IF
+               ELSE         ! these calls are mainly retained to assist in testing
+                  IF ( ln_hpg_ref_ccs ) THEN 
+                     CALL ref_to_tgt_ccs ( jio, jjo, gde3w, rhd, zz_ref, zrhd_ref, zdrhd_k_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) )
+                  ELSE  
+                     CALL ref_to_tgt_cub ( jio, jjo, gde3w, rhd, zz_ref, zrhd_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) )  
+                  END IF
+               END IF
             END IF 
- 
+
             IF ( ln_dbg_hpg ) CALL dbg_3dr( '3. z_rhd_pmr', z_rhd_pmr(:,:,:,jr) ) 
 
@@ -2045,5 +2065,8 @@
       REAL(wp), DIMENSION(A2D(nn_hls),jpk),   INTENT(OUT) :: p_fld_tgt_ref ! target minus reference on target grid         
 
+!---------------------------------------------------------------------------------------------------------------
+
       INTEGER,  DIMENSION(A2D(nn_hls),jpk) :: jk_ref_for_tgt   ! reference index for interpolation to target grid; target lies between jk_ref and jk_ref-1.
+      LOGICAL,  DIMENSION(A2D(nn_hls),jpk) :: ll_tgt_off_cen   ! T => off-centred interpolation (values not used in this routine) 
 
 !---------------------------------------------------------------------------------------------------------------
@@ -2068,5 +2091,5 @@
 
 ! find jk_ref_for_tgt (bounding levels on reference grid for each target point
-      CALL loc_ref_tgt ( ll_ccs, ki_off_tgt, kj_off_tgt, p_dep_tgt, p_z_ref, kk_bot_ref, jk_ref_for_tgt )
+      CALL loc_ref_tgt ( ll_ccs, ki_off_tgt, kj_off_tgt, p_dep_tgt, p_z_ref, kk_bot_ref, jk_ref_for_tgt, ll_tgt_off_cen )
 
       DO_3D( 0, 0, 0, 0, 1, jpk-1 ) 
@@ -2106,4 +2129,88 @@
 
    END SUBROUTINE ref_to_tgt_cub
+
+!---------------------------------------------------------------------------------------------------------------
+
+   SUBROUTINE ref_to_tgt_cub_str ( ki_off_tgt, kj_off_tgt, p_dep_tgt, p_fld_tgt, p_z_ref, p_fld_ref, kk_bot_ref, p_fld_tgt_ref) 
+       
+      INTEGER,                                INTENT(in)  :: ki_off_tgt    ! offset of points in target array in i-direction   
+      INTEGER,                                INTENT(in)  :: kj_off_tgt    ! offset of points in target array in j-direction   
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk),   INTENT(in)  :: p_dep_tgt     ! depths of target profiles       
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk),   INTENT(in)  :: p_fld_tgt     ! field values on the target grid 
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk),   INTENT(in)  :: p_z_ref       ! heights of reference  profiles       
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk),   INTENT(in)  :: p_fld_ref     ! field values to be interpolated (in the vertical) on reference grid
+      INTEGER,  DIMENSION(A2D(nn_hls)),       INTENT(in)  :: kk_bot_ref    ! bottom levels in the reference profile
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk),   INTENT(OUT) :: p_fld_tgt_ref ! target minus reference on target grid         
+
+!---------------------------------------------------------------------------------------------------------------
+
+      INTEGER,  DIMENSION(A2D(nn_hls),jpk) :: jk_ref_for_tgt   ! reference index for interpolation to target grid; target lies between jk_ref and jk_ref-1.
+      LOGICAL,  DIMENSION(A2D(nn_hls),jpk) :: ll_tgt_off_cen   ! T => off-centred interpolation (values not used in this routine) 
+
+!---------------------------------------------------------------------------------------------------------------
+
+      INTEGER  :: ji, jj, jk                   ! loop indices 
+      INTEGER  :: jkr
+
+      REAL(wp) :: zf_p, zf_0, zf_m, zf_b   ! plus, zero, minus, below ; note that zeta increases with depth ; zero is zeta = 1/2
+      REAL(wp) :: zz_p, zz_0, zz_m, zz_b  
+      REAL(wp) :: zs_p, zs_0, zs_m, zs_b
+      REAL(wp) :: zz_p_m, zz_p_b, zz_m_b 
+      REAL(wp) :: za_1, za_2, za_3   
+      
+      REAL(wp) :: zz_tgt_lcl, zz_tgt_sq, zfld_ref_on_tgt 
+      LOGICAL  :: ll_ccs  
+
+!-----------------------------------------------------------------------------------------------------------------
+
+      ll_ccs = .FALSE.      ! set for the call to loc_ref_tgt 
+
+! find jk_ref_for_tgt (bounding levels on reference grid for each target point
+      CALL loc_ref_tgt ( ll_ccs, ki_off_tgt, kj_off_tgt, p_dep_tgt, p_z_ref, kk_bot_ref, jk_ref_for_tgt, ll_tgt_off_cen )
+
+      DO_3D( 0, 0, 0, 0, 1, jpk-1 ) 
+
+! it would probably be better computationally for fld_ref to have the jk index first. 
+
+!!! jkr >= 2 and p_fld_ref has jk = 0 available   
+ 
+         jkr  = jk_ref_for_tgt(ji,jj,jk)	 
+         zf_b = p_fld_ref(ji,jj,jkr-2)
+	 zf_m = p_fld_ref(ji,jj,jkr-1)
+	 zf_0 = p_fld_ref(ji,jj,jkr  )
+	 zf_p = p_fld_ref(ji,jj,jkr+1)
+
+         zz_b = p_z_ref( ji, jj, jkr-2) 
+         zz_m = p_z_ref( ji, jj, jkr-1)
+         zz_0 = p_z_ref( ji, jj, jkr  )   
+         zz_b = zz_b - zz_0 
+         zz_m = zz_m - zz_0  
+         zz_p = p_z_ref( ji, jj, jkr+1) - zz_0
+
+         zz_p_m = zz_p - zz_m 
+         zz_p_b = zz_p - zz_b 
+         zz_m_b = zz_m - zz_b 
+
+         zs_0 = - zf_0 / ( zz_p * zz_m   * zz_b   )    
+
+         zs_p =   zf_p / ( zz_p * zz_p_m * zz_p_b ) 
+         zs_b =   zf_b / ( zz_b * zz_p_b * zz_m_b ) 
+         zs_m = - zf_m / ( zz_m * zz_p_m * zz_m_b ) 
+
+         za_1 =    zz_m*zz_b *zs_p +  zz_p*zz_b *zs_m +  zz_p*zz_m *zs_b + (zz_p*zz_m+zz_m*zz_b+zz_p*zz_b)*zs_0  
+         za_2 = - (zz_m+zz_b)*zs_p - (zz_p+zz_b)*zs_m - (zz_p+zz_m)*zs_b - (zz_p+zz_m+zz_b)*zs_0 
+         za_3 = zs_p + zs_m + zs_0 + zs_b
+
+         zz_tgt_lcl = - p_dep_tgt( ji+ki_off_tgt, jj+kj_off_tgt, jk ) - zz_0
+         zz_tgt_sq = zz_tgt_lcl*zz_tgt_lcl
+	 
+         zfld_ref_on_tgt = zf_0 + za_1*zz_tgt_lcl + za_2*zz_tgt_sq + za_3*zz_tgt_lcl*zz_tgt_sq 
+	 
+! when zfld_ref_on_tgt is commented out in the next line, the results for hpg_djr should agree with those for hpg_djc.   
+         p_fld_tgt_ref(ji, jj, jk) = p_fld_tgt(ji+ki_off_tgt, jj+kj_off_tgt, jk) - zfld_ref_on_tgt
+
+      END_3D   
+
+   END SUBROUTINE ref_to_tgt_cub_str
 
 !-----------------------------------------------------------------------------------------------------------------
@@ -2123,21 +2230,20 @@
 !-----------------------------------------------------------------------------------------------------------------
       INTEGER,  DIMENSION(A2D(nn_hls),jpk) :: jk_ref_for_tgt   ! reference index for interpolation to target grid
+      LOGICAL,  DIMENSION(A2D(nn_hls),jpk) :: ll_tgt_off_cen   ! T => off-centred interpolation (values not used in this routine) 
 
       INTEGER  :: ji, jj, jk                 ! DO loop indices  
       INTEGER  :: jkr
-      REAL(wp) :: z_r_6                      ! constant
       REAL(wp) :: zz_tgt_lcl, zz_ref_jkrm1, zz_ref_jkr, zeta, zetasq
       REAL(wp) :: zd_dif, zd_ave, zf_dif, zf_ave
       REAL(wp) :: zcoef_0, zcoef_1, zcoef_2, zcoef_3    
       REAL(wp) :: zfld_ref_on_tgt
-      LOGICAL  :: ll_ccs                     ! set to .FALSE. for the call to loc_ref_tgt  
+      LOGICAL  :: ll_ccs                     ! set to .TRUE. for the call to loc_ref_tgt  
 
 !-----------------------------------------------------------------------------------------------------------------
 
       ll_ccs = .TRUE.      ! set for the call to loc_ref_tgt 
-      z_r_6  = 1._wp / 6._wp 
 
 ! find jk_ref_for_tgt (bounding levels on reference grid for each target point
-      CALL loc_ref_tgt ( ll_ccs, ki_off_tgt, kj_off_tgt, pdep_tgt, pz_ref, kk_bot_ref, jk_ref_for_tgt )
+      CALL loc_ref_tgt ( ll_ccs, ki_off_tgt, kj_off_tgt, pdep_tgt, pz_ref, kk_bot_ref, jk_ref_for_tgt, ll_tgt_off_cen )
 
       DO_3D( 0, 0, 0, 0, 1, jpk-1 ) 
@@ -2182,5 +2288,261 @@
 !-----------------------------------------------------------------------------------------------------------------
 
-   SUBROUTINE loc_ref_tgt (ll_ccs, ki_off_tgt, kj_off_tgt, p_dep_tgt, p_z_ref, kk_bot_ref, kk_ref_for_tgt ) 
+   SUBROUTINE ref_to_tgt_ccs_str ( ki_off_tgt, kj_off_tgt, pdep_tgt, pfld_tgt, pz_ref, pfld_ref, pdfld_k_ref, kk_bot_ref, pfld_tgt_ref) 
+       
+      INTEGER,                      INTENT(in)  :: ki_off_tgt    ! offset of points in target array in i-direction   
+      INTEGER,                      INTENT(in)  :: kj_off_tgt    ! offset of points in target array in j-direction   
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk),   INTENT(in)  :: pdep_tgt      ! depths of target profiles       
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk),   INTENT(in)  :: pfld_tgt      ! field values on the target grid 
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk),   INTENT(in)  :: pz_ref        ! heights of reference  profiles       
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk),   INTENT(in)  :: pfld_ref      ! reference field values
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk),   INTENT(in)  :: pdfld_k_ref   ! harmonic averages of vertical differences of reference field
+      INTEGER,  DIMENSION(A2D(nn_hls)),       INTENT(in)  :: kk_bot_ref    ! bottom levels in the reference profile
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk),   INTENT(OUT) :: pfld_tgt_ref  ! target minus reference on target grid         
+
+!-----------------------------------------------------------------------------------------------------------------
+      INTEGER,  DIMENSION(A2D(nn_hls),jpk) :: jk_ref_for_tgt   ! reference index for interpolation to target grid
+      LOGICAL,  DIMENSION(A2D(nn_hls),jpk) :: ll_tgt_off_cen   ! T => off-centred interpolation (values not used in this routine) 
+
+      INTEGER  :: ji, jj, jk                 ! DO loop indices  
+      INTEGER  :: jkr
+      REAL(wp) :: zz_tgt_lcl, zz_tgt_sq
+
+      REAL(wp) :: zeta_p, zeta_0, zeta_m, zeta_b   ! plus, zero, minus, below (note that zeta increases with depth) zero is zeta = 1/2
+      REAL(wp) :: zz_p, zz_0, zz_m, zz_b  
+      REAL(wp) :: zs_p, zs_0, zs_m, zs_b
+      REAL(wp) :: zz_p_m, zz_p_b, zz_m_b 
+      REAL(wp) :: za_1, za_2, za_3   
+      REAL(wp) :: zeta, zetasq
+
+      REAL(wp) :: zd_dif, zd_ave, zf_dif, zf_ave
+      REAL(wp) :: zcoef_0, zcoef_1, zcoef_2, zcoef_3    
+      REAL(wp) :: zfld_ref_on_tgt
+      LOGICAL  :: ll_ccs                     ! set to .FALSE. for the call to loc_ref_tgt  
+
+!-----------------------------------------------------------------------------------------------------------------
+
+      ll_ccs = .TRUE.      ! set for the call to loc_ref_tgt 
+
+! find jk_ref_for_tgt (bounding levels on reference grid for each target point
+      CALL loc_ref_tgt ( ll_ccs, ki_off_tgt, kj_off_tgt, pdep_tgt, pz_ref, kk_bot_ref, jk_ref_for_tgt, ll_tgt_off_cen )
+
+      zeta_p =  1.5_wp
+      zeta_0 =  0.5_wp
+      zeta_m = -0.5_wp
+      zeta_b = -1.5_wp
+
+      DO_3D( 0, 0, 0, 0, 1, jpk-1 ) 
+
+         jkr = jk_ref_for_tgt( ji, jj, jk )
+
+         zz_0 = pz_ref( ji, jj, jkr  )   
+         zz_m = pz_ref( ji, jj, jkr-1) - zz_0
+
+         IF ( jkr .NE. 2) THEN 
+            zz_b = pz_ref( ji, jj, jkr-2) - zz_0  
+         ELSE 
+            zz_b = 2._wp * zz_m       
+         END IF 
+
+         IF ( ll_tgt_off_cen(ji,jj,jk) ) THEN 
+            zz_p = - zz_m            
+	 ELSE  
+            zz_p = pz_ref( ji, jj, jkr+1) - zz_0
+         END IF 
+
+         zz_p_m = zz_p - zz_m 
+         zz_p_b = zz_p - zz_b 
+         zz_m_b = zz_m - zz_b 
+
+         zs_0 = - zeta_0 / ( zz_p * zz_m   * zz_b   )    
+
+         zs_p =   zeta_p / ( zz_p * zz_p_m * zz_p_b ) 
+         zs_m = - zeta_m / ( zz_m * zz_p_m * zz_m_b ) 
+         zs_b =   zeta_b / ( zz_b * zz_p_b * zz_m_b ) 
+
+         za_1 =     zz_m*zz_b*zs_p +   zz_p*zz_b*zs_m +   zz_p*zz_m*zs_b + (zz_p*zz_m+zz_m*zz_b+zz_p*zz_b)*zs_0  
+         za_2 = - (zz_m+zz_b)*zs_p - (zz_p+zz_b)*zs_m - (zz_p+zz_m)*zs_b - (zz_p+zz_m+zz_b)*zs_0 
+         za_3 = zs_p + zs_m + zs_0 + zs_b
+
+         zz_tgt_lcl =  - pdep_tgt( ji+ki_off_tgt, jj+kj_off_tgt, jk ) - zz_0
+         zz_tgt_sq = zz_tgt_lcl*zz_tgt_lcl
+	 
+         zeta = zeta_0 + za_1*zz_tgt_lcl + za_2*zz_tgt_sq + za_3*zz_tgt_lcl*zz_tgt_sq 
+         zetasq = zeta*zeta
+
+         zd_dif =            pdfld_k_ref(ji,jj,jkr) - pdfld_k_ref(ji,jj,jkr-1)  
+         zd_ave = 0.5_wp * ( pdfld_k_ref(ji,jj,jkr) + pdfld_k_ref(ji,jj,jkr-1) )   
+
+         zf_dif =            pfld_ref(ji,jj,jkr) - pfld_ref(ji,jj,jkr-1)  
+         zf_ave = 0.5_wp * ( pfld_ref(ji,jj,jkr) + pfld_ref(ji,jj,jkr-1) )   
+
+         zcoef_0 =          zf_ave - 0.125_wp * zd_dif
+         zcoef_1 = 1.5_wp * zf_dif - 0.5_wp   * zd_ave
+         zcoef_2 = 0.5_wp * zd_dif 
+         zcoef_3 = 2.0_wp * zd_ave - 2._wp    * zf_dif   
+	 
+         zfld_ref_on_tgt = zcoef_0 + zeta*zcoef_1 + zetasq * ( zcoef_2 + zeta*zcoef_3 )  
+
+! when zfld_ref_on_tgt is commented out in the next line, the results for hpg_djr should agree with those for hpg_djc.   
+ 
+         pfld_tgt_ref(ji, jj, jk) = pfld_tgt(ji+ki_off_tgt, jj+kj_off_tgt, jk) - zfld_ref_on_tgt  
+
+!         IF ( ln_dbg_hpg .AND. lwp .AND. ji == ki_dbg_cen .AND. jj == kj_dbg_cen .AND. jk == kk_dbg_cen ) THEN
+!            WRITE(numout,*) ' zeta, zd_dif, zd_ave, zf_dif, zf_ave = ', zeta, zd_dif, zd_ave, zf_dif, zf_ave
+!            WRITE(numout,*) ' zz_tgt_lcl, jkr, zz_ref_jkr, zz_ref_jkrm1 =', zz_tgt_lcl, jkr, zz_ref_jkr, zz_ref_jkrm1 
+!            WRITE(numout,*) ' pfld_ref(ji,jj,jkr), pfld_ref(ji,jj,jkr-1) = ', pfld_ref(ji,jj,jkr), pfld_ref(ji,jj,jkr-1)
+!	    WRITE(numout,*) ' zfld_ref_on_tgt = ', zfld_ref_on_tgt 
+!	    WRITE(numout,*) ' pfld_tgt(ji+ki_off_tgt, jj+kj_off_tgt, jk) = ', pfld_tgt(ji+ki_off_tgt, jj+kj_off_tgt, jk)
+!         END IF             
+	   
+      END_3D   
+
+      IF ( ln_dbg_hpg ) CALL dbg_3dr( 'ref_to_tgt_ccs_str: pfld_tgt_ref', pfld_tgt_ref ) 
+
+   END SUBROUTINE ref_to_tgt_ccs_str
+
+!-----------------------------------------------------------------------------------------------------------------
+
+   SUBROUTINE ref_to_tgt_ccs_str_off ( ki_off_tgt, kj_off_tgt, pdep_tgt, pfld_tgt, pz_ref, pfld_ref, pdfld_k_ref, kk_bot_ref, pfld_tgt_ref) 
+
+! Coded for a stretched grid and to perform off-centred pure cubic interpolation at the upper and lower boundaries 
+       
+      INTEGER,                      INTENT(in)  :: ki_off_tgt    ! offset of points in target array in i-direction   
+      INTEGER,                      INTENT(in)  :: kj_off_tgt    ! offset of points in target array in j-direction   
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk),   INTENT(in)  :: pdep_tgt      ! depths of target profiles       
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk),   INTENT(in)  :: pfld_tgt      ! field values on the target grid 
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk),   INTENT(in)  :: pz_ref        ! heights of reference  profiles       
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk),   INTENT(in)  :: pfld_ref      ! reference field values
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk),   INTENT(in)  :: pdfld_k_ref   ! harmonic averages of vertical differences of reference field
+      INTEGER,  DIMENSION(A2D(nn_hls)),       INTENT(in)  :: kk_bot_ref    ! bottom levels in the reference profile
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk),   INTENT(OUT) :: pfld_tgt_ref  ! target minus reference on target grid         
+
+!-----------------------------------------------------------------------------------------------------------------
+      INTEGER,  DIMENSION(A2D(nn_hls),jpk) :: jk_ref_for_tgt   ! reference index for interpolation to target grid
+      LOGICAL,  DIMENSION(A2D(nn_hls),jpk) :: ll_tgt_off_cen   ! T => off-centred interpolation
+
+      INTEGER  :: ji, jj, jk                 ! DO loop indices  
+      INTEGER  :: jkr
+      REAL(wp) :: zz_tgt_lcl, zz_tgt_sq
+
+      REAL(wp) :: zz_p, zz_0, zz_m, zz_b  
+      REAL(wp) :: zeta_p, zeta_0, zeta_m, zeta_b
+      REAL(wp) :: zf_p, zf_0, zf_m, zf_b
+      REAL(wp) :: zs_p, zs_0, zs_m, zs_b
+      REAL(wp) :: zz_p_m, zz_p_b, zz_m_b 
+      REAL(wp) :: za_1, za_2,za_3   
+      REAL(wp) :: zeta, zetasq
+
+      REAL(wp) :: zd_dif, zd_ave, zf_dif, zf_ave
+      REAL(wp) :: zcoef_0, zcoef_1, zcoef_2, zcoef_3    
+      REAL(wp) :: zfld_ref_on_tgt
+      LOGICAL  :: ll_ccs                     ! set to .FALSE. for the call to loc_ref_tgt 
+
+!-----------------------------------------------------------------------------------------------------------------
+
+      ll_ccs = .FALSE.      ! set for the call to loc_ref_tgt
+      
+! find jk_ref_for_tgt (bounding levels on reference grid for each target point) and ll_tgt_off_cen
+      CALL loc_ref_tgt ( ll_ccs, ki_off_tgt, kj_off_tgt, pdep_tgt, pz_ref, kk_bot_ref, jk_ref_for_tgt, ll_tgt_off_cen )
+
+      zeta_p =  1.5_wp
+      zeta_0 =  0.5_wp
+      zeta_m = -0.5_wp
+      zeta_b = -1.5_wp
+
+      DO_3D( 0, 0, 0, 0, 1, jpk-1 )  
+
+         jkr = jk_ref_for_tgt( ji, jj, jk )
+
+         zz_b = pz_ref( ji, jj, jkr-2) 
+         zz_m = pz_ref( ji, jj, jkr-1)
+         zz_0 = pz_ref( ji, jj, jkr  )   
+         zz_b = zz_b - zz_0 
+         zz_m = zz_m - zz_0  
+         zz_p = pz_ref( ji, jj, jkr+1) - zz_0
+
+         zz_tgt_lcl =  - pdep_tgt( ji+ki_off_tgt, jj+kj_off_tgt, jk ) - zz_0
+         zz_tgt_sq = zz_tgt_lcl*zz_tgt_lcl
+
+         IF ( ll_tgt_off_cen(ji,jj,jk) ) THEN 
+
+! off centred pure cubic fit 
+
+            zf_b = pfld_ref(ji,jj,jkr-2)
+            zf_m = pfld_ref(ji,jj,jkr-1)
+            zf_0 = pfld_ref(ji,jj,jkr  )
+            zf_p = pfld_ref(ji,jj,jkr+1)
+
+            zz_p_m = zz_p - zz_m 
+            zz_p_b = zz_p - zz_b 
+            zz_m_b = zz_m - zz_b 
+
+            zs_0 = - zf_0 / ( zz_p * zz_m   * zz_b   )    
+
+            zs_p =   zf_p / ( zz_p * zz_p_m * zz_p_b ) 
+            zs_b =   zf_b / ( zz_b * zz_p_b * zz_m_b ) 
+            zs_m = - zf_m / ( zz_m * zz_p_m * zz_m_b ) 
+
+            za_1 =    zz_m*zz_b *zs_p +  zz_p*zz_b *zs_m +  zz_p*zz_m *zs_b + (zz_p*zz_m+zz_m*zz_b+zz_p*zz_b)*zs_0  
+            za_2 = - (zz_m+zz_b)*zs_p - (zz_p+zz_b)*zs_m - (zz_p+zz_m)*zs_b - (zz_p+zz_m+zz_b)*zs_0 
+            za_3 = zs_p + zs_m + zs_0 + zs_b
+
+            zfld_ref_on_tgt = zf_0 + za_1*zz_tgt_lcl + za_2*zz_tgt_sq + za_3*zz_tgt_lcl*zz_tgt_sq 
+
+         ELSE 
+
+            zz_p_m = zz_p - zz_m 
+            zz_p_b = zz_p - zz_b 
+            zz_m_b = zz_m - zz_b 
+
+            zs_0 = - zeta_0 / ( zz_p * zz_m   * zz_b   )    
+
+            zs_p =   zeta_p / ( zz_p * zz_p_m * zz_p_b ) 
+            zs_m = - zeta_m / ( zz_m * zz_p_m * zz_m_b ) 
+            zs_b =   zeta_b / ( zz_b * zz_p_b * zz_m_b ) 
+
+            za_1 =     zz_m*zz_b*zs_p +   zz_p*zz_b*zs_m +   zz_p*zz_m*zs_b + (zz_p*zz_m+zz_m*zz_b+zz_p*zz_b)*zs_0  
+            za_2 = - (zz_m+zz_b)*zs_p - (zz_p+zz_b)*zs_m - (zz_p+zz_m)*zs_b - (zz_p+zz_m+zz_b)*zs_0 
+            za_3 = zs_p + zs_m + zs_0 + zs_b
+	 
+            zeta = zeta_0 + za_1*zz_tgt_lcl + za_2*zz_tgt_sq + za_3*zz_tgt_lcl*zz_tgt_sq 
+            zetasq = zeta*zeta
+
+            zd_dif =            pdfld_k_ref(ji,jj,jkr) - pdfld_k_ref(ji,jj,jkr-1)  
+            zd_ave = 0.5_wp * ( pdfld_k_ref(ji,jj,jkr) + pdfld_k_ref(ji,jj,jkr-1) )   
+
+            zf_dif =            pfld_ref(ji,jj,jkr) - pfld_ref(ji,jj,jkr-1)  
+            zf_ave = 0.5_wp * ( pfld_ref(ji,jj,jkr) + pfld_ref(ji,jj,jkr-1) )   
+
+            zcoef_0 =          zf_ave - 0.125_wp * zd_dif
+            zcoef_1 = 1.5_wp * zf_dif - 0.5_wp   * zd_ave
+            zcoef_2 = 0.5_wp * zd_dif 
+            zcoef_3 = 2.0_wp * zd_ave - 2._wp    * zf_dif   
+	 
+            zfld_ref_on_tgt = zcoef_0 + zeta*zcoef_1 + zetasq * ( zcoef_2 + zeta*zcoef_3 )  
+
+         END IF 
+
+! when zfld_ref_on_tgt is commented out in the next line, the results for hpg_djr should agree with those for hpg_djc.   
+ 
+         pfld_tgt_ref(ji, jj, jk) = pfld_tgt(ji+ki_off_tgt, jj+kj_off_tgt, jk) - zfld_ref_on_tgt  
+	 
+!         IF ( ln_dbg_hpg .AND. lwp .AND. ji == ki_dbg_cen .AND. jj == kj_dbg_cen .AND. jk == kk_dbg_cen ) THEN
+!            WRITE(numout,*) ' zeta, zd_dif, zd_ave, zf_dif, zf_ave = ', zeta, zd_dif, zd_ave, zf_dif, zf_ave
+!            WRITE(numout,*) ' zz_tgt_lcl, jkr, zz_ref_jkr, zz_ref_jkrm1 =', zz_tgt_lcl, jkr, zz_ref_jkr, zz_ref_jkrm1 
+!            WRITE(numout,*) ' pfld_ref(ji,jj,jkr), pfld_ref(ji,jj,jkr-1) = ', pfld_ref(ji,jj,jkr), pfld_ref(ji,jj,jkr-1)
+!	    WRITE(numout,*) ' zfld_ref_on_tgt = ', zfld_ref_on_tgt 
+!	    WRITE(numout,*) ' pfld_tgt(ji+ki_off_tgt, jj+kj_off_tgt, jk) = ', pfld_tgt(ji+ki_off_tgt, jj+kj_off_tgt, jk)
+!         END IF             
+	   
+      END_3D   
+
+      IF ( ln_dbg_hpg ) CALL dbg_3dr( 'ref_to_tgt_ccs_str_off: pfld_tgt_ref', pfld_tgt_ref ) 
+
+   END SUBROUTINE ref_to_tgt_ccs_str_off
+
+!-----------------------------------------------------------------------------------------------------------------
+
+   SUBROUTINE loc_ref_tgt (ll_ccs, ki_off_tgt, kj_off_tgt, p_dep_tgt, p_z_ref, kk_bot_ref, kk_ref_for_tgt, ll_tgt_off_cen ) 
 
       LOGICAL,                    INTENT(in)   :: ll_ccs           ! true => constrained cubic spline; false => pure cubic fit 
@@ -2191,5 +2553,7 @@
       INTEGER,  DIMENSION(A2D(nn_hls)),     INTENT(in)   :: kk_bot_ref       ! bottom levels in the reference profile
       INTEGER,  DIMENSION(A2D(nn_hls),jpk), INTENT(OUT)  :: kk_ref_for_tgt   ! reference index for interpolation to target grid (the lower point)
-                                                                             ! with jkr = kk_ref_for_tgt(ji,jj,jk) the reference levels are jkr-1 and jkr         
+                                                                             ! with jkr = kk_ref_for_tgt(ji,jj,jk) the reference levels are jkr-1 and jkr
+      LOGICAL,  DIMENSION(A2D(nn_hls),jpk), INTENT(OUT)  :: ll_tgt_off_cen   ! T => target point is off-centred (only applies when ll_ccs is False) 
+       
 !-----------------------------------------------------------------------------------------------------------------
 
@@ -2211,13 +2575,12 @@
 ! 1. Initialisation for search for points on reference grid bounding points on the target grid
 ! the first point on the target grid is assumed to lie between the first two points on the reference grid 
-      IF ( ll_ccs ) THEN
-          jk_init = 2
-      ELSE 
-          jk_init = 3      ! for simple cubic use off-centred interpolation near the surface 
-      END IF 
+
+      jk_init = 2  !   for all cases; enforcement of off-centred interpolation is now done at the end of the routine 
 
       jk_ref(:,:) = jk_init
       kk_ref_for_tgt(:,:,1) = jk_init
       jk_tgt(:,:) = 2 
+
+      ll_tgt_off_cen(:,:,:) = .FALSE.     
 
 ! 2. Find kk_ref_for_tgt
@@ -2322,8 +2685,26 @@
 ! This assumes that every sea point has at least 4 levels.   
 
-      IF ( .NOT. ll_ccs ) THEN
+      IF ( ll_ccs ) THEN
          DO_3D(0, 0, 0, 0, 2, jpk-1) 
-            IF ( kk_ref_for_tgt(ji,jj, jk) == kk_bot_ref(ji,jj) )  kk_ref_for_tgt(ji,jj, jk) = kk_ref_for_tgt(ji,jj, jk) - 1
+            IF ( kk_ref_for_tgt(ji,jj,jk) == kk_bot_ref(ji,jj) ) THEN 
+               ll_tgt_off_cen(ji,jj,jk) = .TRUE.            
+            END IF 
          END_3D
+      ELSE 
+         DO_3D(0, 0, 0, 0, 1, jpk-1) 
+            IF ( kk_ref_for_tgt(ji,jj,jk) == 2 ) THEN 
+               ll_tgt_off_cen(ji,jj,jk) = .TRUE.            
+               kk_ref_for_tgt(ji,jj,jk) = 3
+            END IF 
+            IF ( kk_ref_for_tgt(ji,jj,jk) == kk_bot_ref(ji,jj) ) THEN 
+               ll_tgt_off_cen(ji,jj,jk) = .TRUE.            
+               kk_ref_for_tgt(ji,jj,jk) =  kk_bot_ref(ji,jj) - 1
+            END IF 
+         END_3D
+      END IF 
+
+      IF ( lwp .AND. ln_dbg_hpg ) THEN 
+         WRITE( numout,*) 'loc_ref_tgt at end of section 4 ', ki_off_tgt, kj_off_tgt
+         CALL dbg_3di( 'kk_ref_for_tgt', kk_ref_for_tgt ) 
       END IF 
 
@@ -2411,7 +2792,7 @@
 ! 1. Calculate the reference profile from all points in the stencil
 
-         IF ( ln_hpg_ffr_ref ) THEN 
+         IF ( ln_hpg_ref ) THEN 
             CALL calc_rhd_ref(j_uv, jn_hor_pts, gdept(:,:,:,Kmm), zrhd_ref, zz_ref, jk_bot_ref)    
-            IF ( ln_hpg_ffr_ref_ccs ) CALL calc_drhd_k(zrhd_ref, jk_bot_ref, zdrhd_k_ref) 
+            IF ( ln_hpg_ref_ccs ) CALL calc_drhd_k(zrhd_ref, jk_bot_ref, zdrhd_k_ref) 
          END IF 
 
@@ -2434,11 +2815,22 @@
            END IF 
 
-           IF ( ln_hpg_ffr_ref ) THEN     
-	   
-             IF ( ln_hpg_ffr_ref_ccs ) THEN 
-               CALL ref_to_tgt_ccs ( ji_ro, jj_ro, gdept, rhd, zz_ref, zrhd_ref, zdrhd_k_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) )
-             ELSE 
-               CALL ref_to_tgt_cub ( ji_ro, jj_ro, gdept, rhd, zz_ref, zrhd_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) )  
-             END IF 
+           IF ( ln_hpg_ref ) THEN 
+              IF ( ln_hpg_ref_str ) THEN 
+                 IF ( ln_hpg_ref_ccs ) THEN 
+                    IF ( ln_hpg_ref_off ) THEN 
+                       CALL ref_to_tgt_ccs_str_off ( ji_ro, jj_ro, gdept, rhd, zz_ref, zrhd_ref, zdrhd_k_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) )
+                    ELSE
+                       CALL ref_to_tgt_ccs_str ( ji_ro, jj_ro, gdept, rhd, zz_ref, zrhd_ref, zdrhd_k_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) )
+                    END IF 
+                 ELSE  
+                    CALL ref_to_tgt_cub_str ( ji_ro, jj_ro, gdept, rhd, zz_ref, zrhd_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) )  
+                 END IF
+              ELSE       ! these calls are mainly retained to assist in testing 
+                 IF ( ln_hpg_ref_ccs ) THEN 
+                    CALL ref_to_tgt_ccs ( ji_ro, jj_ro, gdept, rhd, zz_ref, zrhd_ref, zdrhd_k_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) )
+                 ELSE  
+                    CALL ref_to_tgt_cub ( ji_ro, jj_ro, gdept, rhd, zz_ref, zrhd_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) )  
+                 END IF
+              END IF
 
            ELSE !  no subtraction of reference profile 
@@ -2448,5 +2840,5 @@
              END_3D
 
-           END IF ! ln_hpg_ffr_ref ) THEN 
+           END IF ! ln_hpg_ref 
 
            DO_3D (0,0,0,0,1,jpk)
