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 15726)
+++ NEMO/branches/UKMO/NEMO_4.0-TRUNK_r14960_HPG/src/OCE/DYN/dynhpg.F90	(revision 15727)
@@ -1,3 +1,6 @@
 MODULE dynhpg
+
+   !!  v2:  takes into account the stretching of the vertical grid in hpg_ffr
+
    !!======================================================================
    !!                       ***  MODULE  dynhpg  ***
@@ -1738,5 +1741,5 @@
 
          IF ( ln_dbg_hpg ) CALL dbg_2dr( '2.1 ht_0', ht_0 ) 
-         CALL calc_rhd_ref(j_uv, jn_hor_pts, zrhd_ref, zz_ref, jk_bot_ref)   ! Uses rhd (IN) to calculate all other fields (OUT)
+         CALL calc_rhd_ref(j_uv, jn_hor_pts, gde3w, zrhd_ref, zz_ref, jk_bot_ref)   ! Uses rhd (IN) to calculate all other fields (OUT)
 
          IF ( ln_dbg_hpg ) THEN 
@@ -2056,12 +2059,14 @@
       
       REAL(wp) :: zz_tgt_lcl, zfld_ref_on_tgt 
+      LOGICAL  :: ll_ccs  
 
 !-----------------------------------------------------------------------------------------------------------------
 
+      ll_ccs = .FALSE.      ! set for the call to loc_ref_tgt 
       z_r_6  = 1._wp / 6._wp 
       z_r_24 = 1._wp / 24._wp 
 
 ! find jk_ref_for_tgt (bounding levels on reference grid for each target point
-      CALL loc_ref_tgt ( 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 )
 
       DO_3D( 0, 0, 0, 0, 1, jpk-1 ) 
@@ -2099,6 +2104,4 @@
 
       END_3D   
-
-      RETURN 	     
 
    END SUBROUTINE ref_to_tgt_cub
@@ -2128,10 +2131,13 @@
       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 
       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 ( 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 )
 
       DO_3D( 0, 0, 0, 0, 1, jpk-1 ) 
@@ -2172,12 +2178,11 @@
       IF ( ln_dbg_hpg ) CALL dbg_3dr( 'ref_to_tgt_ccs: pfld_tgt_ref', pfld_tgt_ref ) 
 
-      RETURN 	     
-
    END SUBROUTINE ref_to_tgt_ccs
 
 !-----------------------------------------------------------------------------------------------------------------
 
-   SUBROUTINE loc_ref_tgt ( ki_off_tgt, kj_off_tgt, p_dep_tgt, p_z_ref, kk_bot_ref, kk_ref_for_tgt ) 
-
+   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 ) 
+
+      LOGICAL,                    INTENT(in)   :: ll_ccs           ! true => constrained cubic spline; false => pure cubic fit 
       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   
@@ -2206,5 +2211,5 @@
 ! 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 ( ln_hpg_djr_ref_ccs ) THEN
+      IF ( ll_ccs ) THEN
           jk_init = 2
       ELSE 
@@ -2288,6 +2293,6 @@
             IF ( ( z_above < z_tgt .AND. jkr > jk_init )  .OR. z_below > (z_tgt + tol_wp) ) THEN 
                IF ( jcount < 10 ) THEN 
-                  WRITE(numout,*) 'loc_ref_tgt: ji, jj, jk, ki_off_tgt, kj_off_tgt, z_tgt, z_below, z_above '  
-                  WRITE(numout,*) ji, jj, jk, ki_off_tgt, kj_off_tgt, z_tgt, z_below, z_above  
+                  WRITE(numout,*) 'loc_ref_tgt: ji, jj, jk, ki_off_tgt, kj_off_tgt, jkr, z_tgt, z_below, z_above '  
+                  WRITE(numout,*) ji, jj, jk, ki_off_tgt, kj_off_tgt, jkr, z_tgt, z_below, z_above  
                END IF 
                jcount = jcount + 1  
@@ -2317,11 +2322,9 @@
 ! This assumes that every sea point has at least 4 levels.   
 
-      IF ( .NOT. ln_hpg_djr_ref_ccs ) THEN
+      IF ( .NOT. 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
          END_3D
       END IF 
-
-      RETURN 
 
    END SUBROUTINE loc_ref_tgt
@@ -2350,6 +2353,8 @@
       
       ! The following fields could probably be made single level or at most 2 level fields 
-      REAL(wp), DIMENSION(A2D(nn_hls),jpk,4) :: z_e3t_pmr, z_depw_pmr  ! corresponding e3t and gdepw pmr profiles 
-      REAL(wp), DIMENSION(A2D(nn_hls),jpk)   :: z_rhd_mid, z_e3t_mid   ! profiles and e3t interpolated to middle of cell (using ccs)  
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk,4) :: z_dept_pmr, z_depw_pmr ! corresponding gdept and gdepw pmr profiles 
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk)   :: z_rhd_mid              ! rhd profiles interpolated to middle of cell (using ccs or cub)  
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk)   :: z_dept_mid, z_depw_mid ! dept and depw similarly interpolated to middle of cell (using ccs or cub)  
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk)   :: z_e3t_mid              ! e3t profiles interpolated to middle of cell; only printed as a diagnostic  
       REAL(wp), DIMENSION(A2D(nn_hls),jpk,3) :: z_ddepw_ij             ! constrained spline horizontal differences in gdepw 
       REAL(wp), DIMENSION(A2D(nn_hls),jpk,3) :: z_p_pmr, z_F_pmr       ! pressures and forces on vertical faces  
@@ -2365,5 +2370,6 @@
       INTEGER  ::   ji_ro, jj_ro, jr_offset     ! offsets for the reference profile 
       INTEGER  ::   jn_hor_pts                  ! number of profiles required in horizontal (4 if cubic interpolarion or 2 if not) 
-
+      INTEGER  ::   j_t_levs, j_w_levs          ! indicators that field passed in is valid on t-levels or w-levels
+      
       REAL(wp) ::   z_r_6                       ! 1/6 
       REAL(wp) ::   zr_a, zr_b, zr_c            ! rhd evaluated at i, i+1/2 and i+1  
@@ -2379,4 +2385,7 @@
       
       z_r_6 = 1.0_wp/6.0_wp
+      
+      j_t_levs = 1      ! this is not very neat - it is duplicated in the calling routines
+      j_w_levs = 2 
 
       IF ( ln_hpg_ffr_hor_ccs .OR. ln_hpg_ffr_hor_cub ) THEN  
@@ -2400,8 +2409,8 @@
          END IF 
 
-! 1. Calculate the referene profile from all points in the stencil
+! 1. Calculate the reference profile from all points in the stencil
 
          IF ( ln_hpg_ffr_ref ) THEN 
-            CALL calc_rhd_ref(j_uv, jn_hor_pts, zrhd_ref, zz_ref, jk_bot_ref)   ! Uses rhd (IN) to calculate all other fields (OUT)  
+            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) 
          END IF 
@@ -2425,15 +2434,15 @@
            END IF 
 
-           IF ( ln_hpg_ffr_ref ) THEN       !!!! MJB  should gde3w below now be gdept ?? !!!!  
-
+           IF ( ln_hpg_ffr_ref ) THEN     
+	   
              IF ( ln_hpg_ffr_ref_ccs ) THEN 
-               CALL ref_to_tgt_ccs ( ji_ro, jj_ro, gde3w, rhd, zz_ref, zrhd_ref, zdrhd_k_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) )
+               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, gde3w, rhd, zz_ref, zrhd_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) )  
+               CALL ref_to_tgt_cub ( ji_ro, jj_ro, gdept, rhd, zz_ref, zrhd_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) )  
              END IF 
 
            ELSE !  no subtraction of reference profile 
 
-             DO_3D (0,0,0,0,1,jpk-1)
+             DO_3D (0,0,0,0,1,jpk)
                z_rhd_pmr(ji,jj,jk,jr) = rhd(ji+ji_ro,jj+jj_ro,jk)
              END_3D
@@ -2442,5 +2451,5 @@
 
            DO_3D (0,0,0,0,1,jpk)
-               z_e3t_pmr (ji,jj,jk,jr) =   e3t(ji+ji_ro,jj+jj_ro,jk,Kmm)
+               z_dept_pmr(ji,jj,jk,jr) = gdept(ji+ji_ro,jj+jj_ro,jk,Kmm)
                z_depw_pmr(ji,jj,jk,jr) = gdepw(ji+ji_ro,jj+jj_ro,jk,Kmm)
            END_3D
@@ -2458,30 +2467,33 @@
 
             IF ( ln_hpg_ffr_hor_ccs ) THEN
-	       CALL calc_mid_ccs(j_uv, aco_bc_rhd_hor, bco_bc_rhd_hor, z_rhd_pmr, z_rhd_mid)  
-	       CALL calc_mid_ccs(j_uv, aco_bc_z_hor,   bco_bc_z_hor,   z_e3t_pmr, z_e3t_mid)  
+	       CALL calc_mid_ccs(j_uv, j_t_levs, aco_bc_rhd_hor, bco_bc_rhd_hor, z_rhd_pmr,  z_rhd_mid)  
+	       CALL calc_mid_ccs(j_uv, j_t_levs, aco_bc_z_hor,   bco_bc_z_hor,   z_dept_pmr, z_dept_mid)  
+	       CALL calc_mid_ccs(j_uv, j_w_levs, aco_bc_z_hor,   bco_bc_z_hor,   z_depw_pmr, z_depw_mid)  
                CALL calc_dz_dij_ccs( j_uv, z_depw_pmr, z_ddepw_ij )
             ELSE ! ln_hpg_ffr_hor_cub
-	       CALL calc_mid_cub(j_uv, aco_bc_rhd_hor, bco_bc_rhd_hor, z_rhd_pmr, z_rhd_mid)  
-	       CALL calc_mid_cub(j_uv, aco_bc_z_hor,   bco_bc_z_hor,   z_e3t_pmr, z_e3t_mid)  
+	       CALL calc_mid_cub(j_uv, j_t_levs, aco_bc_rhd_hor, bco_bc_rhd_hor, z_rhd_pmr,  z_rhd_mid)  
+	       CALL calc_mid_cub(j_uv, j_t_levs, aco_bc_z_hor,   bco_bc_z_hor,   z_dept_pmr, z_dept_mid)  
+	       CALL calc_mid_cub(j_uv, j_w_levs, aco_bc_z_hor,   bco_bc_z_hor,   z_depw_pmr, z_depw_mid)  
                CALL calc_dz_dij_cub( j_uv, z_depw_pmr, z_ddepw_ij )
             END IF 
 
-            DO_3D (0,0,0,0,1,jpk-1)  
+            DO_3D (0,0,0,0,1,jpk)  
                z_rhd_pmr(ji,jj,jk,1) = z_rhd_pmr(ji,jj,jk,2) 
                z_rhd_pmr(ji,jj,jk,2) = z_rhd_mid(ji,jj,jk)
-               z_e3t_pmr(ji,jj,jk,1) = z_e3t_pmr(ji,jj,jk,2) 
-               z_e3t_pmr(ji,jj,jk,2) = z_e3t_mid(ji,jj,jk)
+               z_dept_pmr(ji,jj,jk,1) = z_dept_pmr(ji,jj,jk,2) 
+               z_dept_pmr(ji,jj,jk,2) = z_dept_mid(ji,jj,jk)
+               z_depw_pmr(ji,jj,jk,1) = z_depw_pmr(ji,jj,jk,2) 
+               z_depw_pmr(ji,jj,jk,2) = z_depw_mid(ji,jj,jk)
             END_3D
 
          ELSE !  simple linear interpolation 
 
-            DO_3D (0,0,0,0,1,jpk-1)  
-               z_rhd_pmr(ji,jj,jk,3) = z_rhd_pmr(ji,jj,jk,2) 
-               z_rhd_pmr(ji,jj,jk,2) = 0.5_wp*( z_rhd_pmr(ji,jj,jk,1) + z_rhd_pmr(ji,jj,jk,2) )
-            END_3D
-            DO_3D (0,0,0,0,1,jpk-1)  
-               z_e3t_pmr(ji,jj,jk,1) =          e3t(ji, jj, jk, Kmm)    
-               z_e3t_pmr(ji,jj,jk,2) = 0.5_wp*( e3t(ji, jj, jk, Kmm) + e3t(ji+jio, jj+jjo, jk, Kmm)  )   
-               z_e3t_pmr(ji,jj,jk,3) =                                 e3t(ji+jio, jj+jjo, jk, Kmm)    
+            DO_3D (0,0,0,0,1,jpk)  
+               z_rhd_pmr(ji,jj,jk,3) =  z_rhd_pmr(ji,jj,jk,2) 
+               z_rhd_pmr(ji,jj,jk,2) =  0.5_wp*( z_rhd_pmr(ji,jj,jk,1)  + z_rhd_pmr(ji,jj,jk,2) )
+               z_dept_pmr(ji,jj,jk,3) = z_dept_pmr(ji,jj,jk,2)    
+               z_dept_pmr(ji,jj,jk,2) = 0.5_wp*( z_dept_pmr(ji,jj,jk,1) + z_dept_pmr(ji,jj,jk,2)  )   
+               z_depw_pmr(ji,jj,jk,3) = z_depw_pmr(ji,jj,jk,2)    
+               z_depw_pmr(ji,jj,jk,2) = 0.5_wp*( z_depw_pmr(ji,jj,jk,1) + z_depw_pmr(ji,jj,jk,2)  )   
             END_3D
 	 END IF 
@@ -2494,8 +2506,16 @@
 
          IF ( ln_dbg_hpg ) THEN 
-            CALL dbg_3dr( '3. e3t ',         e3t ) 
             CALL dbg_3dr( '3. z_rhd_pmr: 1', z_rhd_pmr(:,:,:,1) ) 
             CALL dbg_3dr( '3. z_rhd_pmr: 2', z_rhd_pmr(:,:,:,2) ) 
             CALL dbg_3dr( '3. z_rhd_pmr: 3', z_rhd_pmr(:,:,:,3) ) 
+            CALL dbg_3dr( '3. z_depw_pmr: 1', z_depw_pmr(:,:,:,1) ) 
+            CALL dbg_3dr( '3. z_depw_pmr: 2', z_depw_pmr(:,:,:,2) ) 
+            CALL dbg_3dr( '3. z_depw_pmr: 3', z_depw_pmr(:,:,:,3) ) 
+            DO jr = 1, 3
+               DO_3D (0,0,0,0,1,jpk-1)
+                  z_e3t_mid(ji,jj,jk) = z_depw_pmr(ji,jj,jk,jr) - z_depw_pmr(ji,jj,jk+1,jr)
+               END_3D
+               CALL dbg_3dr( '3. z_e3t_mid jr : ', z_e3t_mid(:,:,:) ) 
+            END DO 
          END IF 
 
@@ -2504,5 +2524,5 @@
          IF ( ln_hpg_ffr_vrt_quad ) THEN 
             DO jr = 1, 3 
-               CALL vrt_int_quad(j_mbkt(:,:,jr), z_rhd_pmr(:,:,:,jr), z_e3t_pmr(:,:,:,jr), z_p_pmr(:,:,:,jr), z_F_pmr(:,:,:,jr)) 
+               CALL vrt_int_quad(j_mbkt(:,:,jr), z_rhd_pmr(:,:,:,jr), z_dept_pmr(:,:,:,jr), z_depw_pmr(:,:,:,jr), z_p_pmr(:,:,:,jr), z_F_pmr(:,:,:,jr)) 
             END DO ! jr
 
@@ -2515,8 +2535,8 @@
                DO jk = 1, jpk - 1
                   DO_2D(0,0,0,0) 
-	             z_p_pmr(ji,jj,jk+1,jr) = z_p_pmr(ji,jj,jk,jr) + grav*z_e3t_pmr(ji,jj,jk,jr)*z_rhd_pmr(ji,jj,jk,jr)
+	             z_p_pmr(ji,jj,jk+1,jr) = z_p_pmr(ji,jj,jk,jr) + grav*(z_depw_pmr(ji,jj,jk+1,jr) - z_depw_pmr(ji,jj,jk,jr)) *z_rhd_pmr(ji,jj,jk,jr)
                   END_2D
                   DO_2D(0,0,0,0) 
-                     z_F_pmr(ji,jj,jk,jr) = 0.5_wp*z_e3t_pmr(ji,jj,jk,jr)*( z_p_pmr(ji,jj,jk,jr) + z_p_pmr(ji,jj,jk+1,jr) ) 
+                     z_F_pmr(ji,jj,jk,jr) = 0.5_wp*(z_depw_pmr(ji,jj,jk+1,jr) - z_depw_pmr(ji,jj,jk,jr))*( z_p_pmr(ji,jj,jk,jr) + z_p_pmr(ji,jj,jk+1,jr) ) 
                   END_2D
                END DO ! jk 
@@ -2531,7 +2551,4 @@
             CALL dbg_3dr( '4. z_F_pmr: 1', z_F_pmr(:,:,:,1) ) 
             CALL dbg_3dr( '4. z_F_pmr: 3', z_F_pmr(:,:,:,3) ) 
-            CALL dbg_3dr( ' z_e3t_pmr: 1', z_e3t_pmr(:,:,:,1) )
-            CALL dbg_3dr( ' z_e3t_pmr: 2', z_e3t_pmr(:,:,:,2) )
-            CALL dbg_3dr( ' z_e3t_pmr: 3', z_e3t_pmr(:,:,:,3) )
          END IF 
 
@@ -2587,7 +2604,15 @@
              CALL iom_put( "u_force_west", z_F_pmr(:,:,:,1) )
              CALL iom_put( "u_force_upper", z_F_upp )
+!             CALL iom_put( "mid_x_press", z_p_pmr(:,:,:,2) ) 
+!	     CALL iom_put( "x_mid_slope", z_ddepw_ij(:,:,:,2) ) 
+!             CALL iom_put( "lhs_x_press", z_p_pmr(:,:,:,1) ) 
+!	     CALL iom_put( "x_lhs_slope", z_ddepw_ij(:,:,:,1) ) 
+!             CALL iom_put( "rhs_x_press", z_p_pmr(:,:,:,3) ) 
+!	     CALL iom_put( "x_rhs_slope", z_ddepw_ij(:,:,:,3) ) 
           ELSE 
              CALL iom_put( "v_force_south", z_F_pmr(:,:,:,1) )
              CALL iom_put( "v_force_upper", z_F_upp )
+!             CALL iom_put( "mid_y_press", z_p_pmr(:,:,:,2) ) 
+!	     CALL iom_put( "y_mid_slope", z_ddepw_ij(:,:,:,2) ) 
           END IF 
 
@@ -2598,5 +2623,5 @@
 !------------------------------------------------------------------------------------------
    
-   SUBROUTINE calc_rhd_ref (k_uv, kn_hor, prhd_ref, pz_ref, kk_bot_ref) 
+   SUBROUTINE calc_rhd_ref (k_uv, kn_hor, pdep, prhd_ref, pz_ref, kk_bot_ref) 
 
       !!---------------------------------------------------------------------
@@ -2612,5 +2637,5 @@
       INTEGER                             , INTENT( in )  ::  k_uv        ! 1 for u-vel; 2 for v_vel
       INTEGER                             , INTENT( in )  ::  kn_hor      ! 4 for cubic; 2 for linear
-
+      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::  pdep       ! depths (either gdept or gde3w)
       REAL(wp), DIMENSION(:,:,:), INTENT(out) ::  prhd_ref   ! densities    of reference profile
       REAL(wp), DIMENSION(:,:,:), INTENT(out) ::  pz_ref     ! heights      of   "         " 
@@ -2667,6 +2692,6 @@
       jir = ji_ref(ji,jj) 
       jjr = jj_ref(ji,jj) 
-      prhd_ref (ji,jj,jk) =     rhd(jir, jjr, jk) 
-      pz_ref   (ji,jj,jk) = - gde3w(jir, jjr, jk) 
+      prhd_ref (ji,jj,jk) =    rhd(jir, jjr, jk) 
+      pz_ref   (ji,jj,jk) = - pdep(jir, jjr, jk) 
    END_3D 
 
@@ -2725,5 +2750,5 @@
 !----------------------------------------------------------------------------
  
-   SUBROUTINE calc_mid_ccs( k_uv, p_aco_bc_fld_hor, p_bco_bc_fld_hor, p_fld_pmr, p_fld_mid ) 
+   SUBROUTINE calc_mid_ccs( k_uv, k_lev_type, p_aco_bc_fld_hor, p_bco_bc_fld_hor, p_fld_pmr, p_fld_mid ) 
 
       !!---------------------------------------------------------------------
@@ -2731,5 +2756,4 @@
       !!
       !! ** Method  :   Use constrained cubic spline to interpolate 4 profiles to their central point  
-      !!                This version can only be used to interpolate fields on tracer levels (not w-levels)
       !!                  
       !! ** Action : - set p_fld_mid 
@@ -2737,4 +2761,5 @@
 
       INTEGER                               , INTENT(in)  :: k_uv             ! 1 for u-vel; 2 for v_vel
+      INTEGER                               , INTENT(in)  :: k_lev_type       ! 1 for t-level data; 2 for w-level data; used with check of land/sea mask
       REAL(wp)                              , INTENT(in)  :: p_aco_bc_fld_hor ! a coeff for horizontal bc (von Neumann or linear extrapolation)
       REAL(wp)                              , INTENT(in)  :: p_bco_bc_fld_hor ! b coeff for horizontal bc (von Neumann or linear extrapolation)
@@ -2745,11 +2770,23 @@
 
       INTEGER  ji, jj, jk    ! standard loop indices
-      INTEGER  ::   j_t_levs ! indicator that field passed is on tracer levels 
-             
+
+      INTEGER  ::   j_t_levs, jk_msk          ! indicators that field passed in is valid on t-levels or w-levels, jk level to use with masks
+       
          j_t_levs =  1
 
          DO jk = 1, jpk  
 
-            CALL calc_dfld_pmr_ij( k_uv, j_t_levs, jk, p_aco_bc_fld_hor, p_bco_bc_fld_hor, p_fld_pmr, zz_dfld_ij )  
+            IF ( k_lev_type == j_t_levs ) THEN 
+	       jk_msk = jk 
+            ELSE 
+               IF ( jk == 1 ) THEN 
+                  jk_msk = 1
+	       ELSE 
+	          jk_msk = jk - 1 
+               END IF 
+            END IF 
+	    
+            CALL calc_dfld_pmr_ij( k_uv, k_lev_type, jk, p_aco_bc_fld_hor, p_bco_bc_fld_hor, p_fld_pmr, zz_dfld_ij )  
+
 
 ! first contribution is simple centred average; p_rhd_pmr has 4 points; 2 and 3 are the central ones 
@@ -2760,5 +2797,5 @@
             IF ( k_uv == 1 ) THEN 
                DO_2D( 0, 0, 0, 0)       
-                  IF ( umask(ji-1, jj, jk) > 0.5 .OR. umask(ji+1, jj, jk) > 0.5 ) THEN
+                  IF ( umask(ji-1, jj, jk_msk) > 0.5 .OR. umask(ji+1, jj, jk_msk) > 0.5 ) THEN
                      p_fld_mid(ji,jj,jk) = p_fld_mid(ji,jj,jk) - 0.125_wp * ( zz_dfld_ij(ji,jj,jk,2) - zz_dfld_ij(ji,jj,jk,1) )
                   END IF 	    
@@ -2766,5 +2803,5 @@
             ELSE !  k_uv == 2 
                DO_2D( 0, 0, 0, 0)       
-                  IF ( vmask(ji, jj-1, jk) > 0.5 .OR. vmask(ji, jj+1, jk) > 0.5 ) THEN
+                  IF ( vmask(ji, jj-1, jk_msk) > 0.5 .OR. vmask(ji, jj+1, jk_msk) > 0.5 ) THEN
                      p_fld_mid(ji,jj,jk) = p_fld_mid(ji,jj,jk) - 0.125_wp * ( zz_dfld_ij(ji,jj,jk,2) - zz_dfld_ij(ji,jj,jk,1) )
                   END IF 	    
@@ -2774,9 +2811,15 @@
 	 END DO  ! jk 
 
+         IF ( ln_dbg_hpg ) THEN 
+            CALL dbg_3dr( 'calc_mid_ccs: zz_dfld_ij: 1', zz_dfld_ij(:,:,:,1) ) 
+            CALL dbg_3dr( 'calc_mid_ccs: zz_dfld_ij: 2', zz_dfld_ij(:,:,:,2) ) 
+         END IF 
+
+
    END SUBROUTINE calc_mid_ccs
 
 !----------------------------------------------------------------------------
  
-   SUBROUTINE calc_mid_cub( k_uv, p_aco_bc_fld_hor, p_bco_bc_fld_hor, p_fld_pmr, p_fld_mid ) 
+   SUBROUTINE calc_mid_cub( k_uv, k_lev_type, p_aco_bc_fld_hor, p_bco_bc_fld_hor, p_fld_pmr, p_fld_mid ) 
 
       !!---------------------------------------------------------------------
@@ -2790,4 +2833,5 @@
 
       INTEGER                               , INTENT(in)  :: k_uv             ! 1 for u-vel; 2 for v_vel
+      INTEGER                               , INTENT(in)  :: k_lev_type       ! 1 for t-level data; 2 for w-level data; used with check of land/sea mask
       REAL(wp)                              , INTENT(in)  :: p_aco_bc_fld_hor ! a coeff for horizontal bc (von Neumann or linear extrapolation)
       REAL(wp)                              , INTENT(in)  :: p_bco_bc_fld_hor ! b coeff for horizontal bc (von Neumann or linear extrapolation)
@@ -2802,7 +2846,21 @@
       REAL(wp) z_cor            ! correction to the central value  
 
+      INTEGER  ::   j_t_levs, jk_msk          ! indicators that field passed in is valid on t-levels or w-levels, jk level to use with masks
+       
+         j_t_levs =  1
+
          z_1_16 = 1._wp / 16._wp   ;    z_9_16 = 9._wp / 16._wp        
 
          DO jk = 1, jpk 
+
+            IF ( k_lev_type == j_t_levs ) THEN 
+	       jk_msk = jk 
+            ELSE 
+               IF ( jk == 1 ) THEN 
+                  jk_msk = 1
+	       ELSE 
+	          jk_msk = jk - 1 
+               END IF 
+            END IF 
 
 ! first contribution is simple centred average plus 1/16 of it; p_rhd_pmr has 4 points; 2 and 3 are the central ones 
@@ -2813,10 +2871,10 @@
             IF ( k_uv == 1 ) THEN 
                DO_2D( 0, 0, 0, 0)       
-                  IF ( umask(ji-1, jj, jk) > 0.5  ) THEN
+                  IF ( umask(ji-1, jj, jk_msk) > 0.5  ) THEN
                      z_cor = p_fld_pmr(ji,jj,jk,1)
                   ELSE 	    
                      z_cor = p_fld_pmr(ji,jj,jk,2)
                   END IF 	    
-                  IF ( umask(ji+1, jj, jk) > 0.5  ) THEN
+                  IF ( umask(ji+1, jj, jk_msk) > 0.5  ) THEN
                      z_cor = z_cor + p_fld_pmr(ji,jj,jk,4)
                   ELSE 	    
@@ -2827,10 +2885,10 @@
             ELSE !  k_uv == 2 
                DO_2D( 0, 0, 0, 0)       
-                  IF ( vmask(ji, jj-1, jk) > 0.5 ) THEN  
+                  IF ( vmask(ji, jj-1, jk_msk) > 0.5 ) THEN  
                      z_cor = p_fld_pmr(ji,jj,jk,1)
                   ELSE 	    
                      z_cor = p_fld_pmr(ji,jj,jk,2)
                   END IF 	    
-                  IF ( vmask(ji, jj+1, jk) > 0.5  ) THEN
+                  IF ( vmask(ji, jj+1, jk_msk) > 0.5  ) THEN
                      z_cor = z_cor + p_fld_pmr(ji,jj,jk,4)
                   ELSE 	    
@@ -2864,9 +2922,9 @@
 
       INTEGER  ji, jj, jk    ! standard loop indices
-      INTEGER  jk_msk                                        ! jk level to use for umask and vmask (we need z on the upper and lower faces)  
+      INTEGER  jk_msk        ! jk level to use for umask and vmask (we need z on the upper and lower faces)  
       INTEGER  j_w_levs      ! indicator for data on w-levels   
-
+      
          j_w_levs = 2 
-               
+                   
          DO jk = 1, jpk 
 
@@ -2878,5 +2936,5 @@
             END_2D 
 
-! set the central element if it is not an isolated velocity cell; this evaluates dz_dij at zeta = 0 ; that is f^{(1)}; and that is given by SMcW(5.8)  
+! set the central element if it is not an isolated velocity cell 
             IF ( jk == 1 ) THEN 
                jk_msk = 1
@@ -3120,5 +3178,5 @@
 !------------------------------------------------------------------------------------------
 
-   SUBROUTINE vrt_int_quad(k_mbkt, p_rhd, p_e3t, p_p, p_F) 
+   SUBROUTINE vrt_int_quad(k_mbkt, p_rhd, p_dept, p_depw, p_p, p_F) 
 
       !!---------------------------------------------------------------------
@@ -3126,5 +3184,5 @@
       !!
       !! ** Method  :   Use quadratic reconstruction of p_rhd (density profile) to calculate pressure profile and 
-      !!                forces on vertical faces between w levels   
+      !!                forces on vertical faces between w levels. Formulated for variable vertical grid spacing.    
       !!                  
       !! ** Action : - set p_p and p_F
@@ -3133,67 +3191,116 @@
       INTEGER,  DIMENSION(A2D(nn_hls)),     INTENT(in)  :: k_mbkt   ! bottom level
       REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in)  :: p_rhd    ! densities on tracer levels
-      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in)  :: p_e3t    ! layer thickness between w levels
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in)  :: p_dept   ! depths of T points 
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in)  :: p_depw   ! depths of w points (top & bottom of T cells)  
       REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(out) :: p_p      ! pressures on w levels
       REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(out) :: p_F      ! force on the vertical face between w levels
 
       INTEGER  ::   ji, jj, jk                  ! dummy loop indices
-      REAL(wp) ::   z_2_3, z_4_3, z_8_3         ! 2/3 and 4/3 and 8/3
-      REAL(wp) ::   z_22_3, z_28_3              ! 22/3 and 28/3 
-
-      REAL(wp) ::   zr_a, zr_b, zr_c            ! rhd evaluated at i, i+1/2 and i+1  
-      REAL(wp) ::   za_0, za_1, za_2            ! polynomial coefficients
-      REAL(wp) ::   ze3t                        ! 
-
-      z_2_3  = 2.0_wp/3.0_wp
-      z_4_3  = 4.0_wp/3.0_wp
-      z_8_3  = 8.0_wp/3.0_wp
-      z_22_3 = 22.0_wp/3.0_wp
-      z_28_3 = 28.0_wp/3.0_wp 
+      REAL(wp) ::   z_1_3                       ! 1/3  
+
+      INTEGER, DIMENSION(A2D(nn_hls))  :: jklev
+      REAL(wp), DIMENSION(A2D(nn_hls)) :: zr_a, zr_0, zr_b   ! rhd values (r) at points a (above), 0 (central), b (below) 
+      REAL(wp), DIMENSION(A2D(nn_hls)) :: zz_a, zz_0, zz_b   ! heights at points a (above), 0 (central), b (below) 
+      REAL(wp), DIMENSION(A2D(nn_hls)) :: zz_upp, zz_low     ! upper & lower w levels 
+      INTEGER  ::   jkt                                 ! level being worked on 
+      REAL(wp) ::   zet_a, zet_0, zet_b                 ! zeta values for above (a), central (0) and below (b) 
+      REAL(wp) ::   zr_d_a, zr_d_0, zr_d_b              ! reciprocals of denominators for terms involving r_a, r_0 and r_b 
+      REAL(wp) ::   za_0, za_1, za_2                    ! quadratic polynomial coefficients for r (rhd)  
+      REAL(wp) ::   zet_upp, zet_low                    ! zeta for upper and lower w levels  
+      REAL(wp) ::   zet_upp_sq, zet_low_sq              ! squares of above values 
+      REAL(wp) ::   zp_upp, zp_low                      ! contributions to the pressure 
+      REAL(wp) ::   zf_uml                              ! force "upper minus lower" 
+      
+      z_1_3  = 1.0_wp/3.0_wp
 
 ! Integrate densities in the vertical using quadratic fits to the densities 
 ! zr_a is r_{-1}  zr_b is r_1 ; zr_c is r_3  
 
+      DO_2D(0,0,0,0) 
+         p_p(ji,jj,1) = 0.0_wp
+      END_2D
+
 ! one-sided quadratic at the upper boundary. pressure is zero at the upper surface
-      DO_2D(0,0,0,0) 
-	 zr_a = p_rhd(ji,jj,1) 
-	 zr_b = p_rhd(ji,jj,2) 
-	 zr_c = p_rhd(ji,jj,3) 
-	 ze3t = p_e3t(ji,jj,1)
-         za_0 = 0.125_wp * ( - zr_c + 6._wp*zr_b + 3._wp*zr_a ) 
-         za_1 = 0.5_wp   * ( zr_b - zr_a )  
-         za_2 = 0.125_wp * (   zr_c - 2._wp*zr_b + zr_a ) 
-         p_p(ji,jj,1) = 0.0_wp	  
-         p_p(ji,jj,2) =        grav*ze3t*     ( za_0 -       za_1 + z_4_3*za_2 )  
-         p_F(ji,jj,1) = 0.5_wp*grav*ze3t*ze3t*( za_0 - z_4_3*za_1 + 2._wp*za_2 ) 
-      END_2D 
-
-! centred quadratic in the interior 
-      DO_3D(0,0,0,0,2,jpk-1) 
-	 zr_a = p_rhd(ji,jj,jk-1) 
-	 zr_b = p_rhd(ji,jj,jk) 
-	 zr_c = p_rhd(ji,jj,jk+1) 
-	 ze3t = p_e3t(ji,jj,jk)
-         za_0 = 0.125_wp * ( - zr_c + 6._wp*zr_b + 3._wp*zr_a ) 
-         za_1 = 0.5_wp   * ( zr_b - zr_a )  
-         za_2 = 0.125_wp * (   zr_c - 2._wp*zr_b + zr_a )  
-         p_p(ji,jj,jk+1) =      p_p(ji,jj,jk) +     grav*ze3t*     ( za_0 +       za_1 + z_4_3*za_2 )  
-         p_F(ji,jj,jk  ) = ze3t*p_p(ji,jj,jk) + 0.5*grav*ze3t*ze3t*( za_0 + z_2_3*za_1 + z_2_3*za_2 ) 
-      END_3D 
-
-! one-sided quadratic at the lower boundary
-      DO_2D(0,0,0,0) 
-         jk   = k_mbkt(ji,jj) 
-	 zr_a = p_rhd(ji,jj,jk-2) 
-	 zr_b = p_rhd(ji,jj,jk-1) 
-	 zr_c = p_rhd(ji,jj,jk) 
-	 ze3t = p_e3t(ji,jj,jk)
-         za_0 = 0.125_wp * ( - zr_c + 6._wp*zr_b + 3._wp*zr_a ) 
-         za_1 = 0.5_wp   * ( zr_b - zr_a )  
-         za_2 = 0.125_wp * (   zr_c - 2._wp*zr_b + zr_a )  
-         p_p(ji,jj,jk+1) =      p_p(ji,jj,jk) +     grav*ze3t*     ( za_0 +    3.*za_1 + z_28_3*za_2 )  
-         p_F(ji,jj,jk  ) = ze3t*p_p(ji,jj,jk) + 0.5*grav*ze3t*ze3t*( za_0 + z_8_3*za_1 + z_22_3*za_2 ) 
-      END_2D 
-
-   RETURN    
+      DO jk = 1, jpk 
+      
+         IF (jk == 1) THEN        !  perform off-centred calculations for the first level 
+            DO_2D(0,0,0,0) 
+               jklev(ji,jj)  = jk 
+	       zr_a(ji,jj)   =   p_rhd(ji,jj,1) 
+	       zr_0(ji,jj)   =   p_rhd(ji,jj,2) 
+	       zr_b(ji,jj)   =   p_rhd(ji,jj,3) 
+	       zz_a(ji,jj)   = - p_dept(ji,jj,1)
+	       zz_0(ji,jj)   = - p_dept(ji,jj,2)
+               zz_b(ji,jj)   = - p_dept(ji,jj,3)
+               zz_upp(ji,jj) = - p_depw(ji,jj,1) 
+	       zz_low(ji,jj) = - p_depw(ji,jj,2) 
+            END_2D
+
+         ELSE IF ( jk < jpk ) THEN 
+            DO_2D(0,0,0,0) 
+               jklev(ji,jj)  = jk 
+	       zr_a(ji,jj)   =   p_rhd(ji,jj,jk-1) 
+	       zr_0(ji,jj)   =   p_rhd(ji,jj,jk) 
+	       zr_b(ji,jj)   =   p_rhd(ji,jj,jk+1) 
+	       zz_a(ji,jj)   = - p_dept(ji,jj,jk-1)
+	       zz_0(ji,jj)   = - p_dept(ji,jj,jk)
+               zz_b(ji,jj)   = - p_dept(ji,jj,jk+1)
+               zz_upp(ji,jj) = - p_depw(ji,jj,jk) 
+	       zz_low(ji,jj) = - p_depw(ji,jj,jk+1) 
+            END_2D
+
+         ELSE IF ( jk == jpk ) THEN    ! repeat calculations at bottom level using off-centred quadratic
+
+            DO_2D(0,0,0,0) 
+               jkt          = k_mbkt(ji,jj) 
+               jklev(ji,jj) = jkt 
+	       zr_a(ji,jj)   =   p_rhd(ji,jj,jkt-2) 
+	       zr_0(ji,jj)   =   p_rhd(ji,jj,jkt-1) 
+	       zr_b(ji,jj)   =   p_rhd(ji,jj,jkt) 
+	       zz_a(ji,jj)   = - p_dept(ji,jj,jkt-2)
+	       zz_0(ji,jj)   = - p_dept(ji,jj,jkt-1)
+               zz_b(ji,jj)   = - p_dept(ji,jj,jkt)
+               zz_upp(ji,jj) = - p_depw(ji,jj,jkt) 
+	       zz_low(ji,jj) = - p_depw(ji,jj,jkt+1) 
+            END_2D
+
+         END IF  ! jk 
+	 
+         DO_2D(0,0,0,0) 
+            jkt  = jklev(ji,jj) 
+
+            zet_a = - ( zz_a(ji,jj) - zz_0(ji,jj) )  
+            zet_0 = 0.0_wp
+	    zet_b = - ( zz_b(ji,jj) - zz_0(ji,jj) )
+
+            zr_d_a = 1.0_wp / ( zet_a*(zet_a - zet_b) ) 
+            zr_d_0 = 1.0_wp / ( zet_a*zet_b ) 
+	    zr_d_b = 1.0_wp / ( zet_b*(zet_b - zet_a) )
+
+            za_0 = zr_0(ji,jj)  
+            za_1 = - zr_d_0*zr_0(ji,jj)*(zet_a+zet_b) - zr_d_a*zr_a(ji,jj)*zet_b - zr_d_b*zr_b(ji,jj)*zet_a  
+            za_2 =   zr_d_0*zr_0(ji,jj)               + zr_d_a*zr_a(ji,jj)       + zr_d_b*zr_b(ji,jj) 
+	    
+            zet_upp = - ( zz_upp(ji,jj) - zz_0(ji,jj) ) 
+            zet_low = - ( zz_low(ji,jj) - zz_0(ji,jj) ) 
+
+            zet_upp_sq = zet_upp*zet_upp   ;  zet_low_sq = zet_low*zet_low
+	    
+	    zp_upp = za_0*zet_upp + 0.5_wp*za_1*zet_upp_sq + z_1_3*za_2*zet_upp_sq*zet_upp  
+            zp_low = za_0*zet_low + 0.5_wp*za_1*zet_low_sq + z_1_3*za_2*zet_low_sq*zet_low 
+
+            p_p(ji,jj,jkt+1) = p_p(ji,jj,jkt) + grav * ( zp_low - zp_upp ) 
+
+            zf_uml =          za_0 * ( 0.5_wp *(zet_low_sq            - zet_upp_sq)            - zet_upp           *(zet_low-zet_upp) ) &
+            &        + 0.5_wp*za_1 * ( z_1_3  *(zet_low_sq*zet_low    - zet_upp_sq*zet_upp)    - zet_upp_sq        *(zet_low-zet_upp) ) &
+            &        + z_1_3 *za_2 * ( 0.25_wp*(zet_low_sq*zet_low_sq - zet_upp_sq*zet_upp_sq) - zet_upp_sq*zet_upp*(zet_low-zet_upp) ) 
+ 
+! the force on the face is positive (the opposite sign from my notes) 
+            p_F(ji,jj,jkt) = p_p(ji,jj,jkt)*( zz_upp(ji,jj) - zz_low(ji,jj) ) + grav * zf_uml   
+
+         END_2D 
+
+      END DO  ! jk  
+
    END SUBROUTINE vrt_int_quad
    
