Index: vendor/nemo/current/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- vendor/nemo/current/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 2)
+++ vendor/nemo/current/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 4)
@@ -14,6 +14,6 @@
    !!             -   !  2005-11  (G. Madec) style & small optimisation
    !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
-   !!            3.4  !  2011-11  (A. Coward, H. Liu) introduction of prj scheme; 
-   !!                 !           suppression of hel, wdj and rot options
+   !!            3.4  !  2011-11  (H. Liu) hpg_prj: Original code for s-coordinates
+   !!                 !           (A. Coward) suppression of hel, wdj and rot options
    !!----------------------------------------------------------------------
 
@@ -61,5 +61,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id: dynhpg.F90 3186 2011-11-27 08:16:19Z smasson $
+   !! $Id: dynhpg.F90 3294 2012-01-28 16:44:18Z rblod $
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -76,6 +76,4 @@
       !!             - Save the trend (l_trddyn=T)
       !!----------------------------------------------------------------------
-      USE oce, ONLY:   tsa                         ! (tsa) used as 2 3D workspace
-      !!
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv
@@ -85,7 +83,5 @@
       !
       IF( l_trddyn ) THEN                    ! Temporary saving of ua and va trends (l_trddyn)
-         ztrdu => tsa(:,:,:,1) 
-         ztrdv => tsa(:,:,:,2) 
-         !
+         CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv )
          ztrdu(:,:,:) = ua(:,:,:)  
          ztrdv(:,:,:) = va(:,:,:) 
@@ -104,4 +100,5 @@
          ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
          CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_hpg, 'DYN', kt )
+         CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv )
       ENDIF          
       !
@@ -191,6 +188,4 @@
       !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
       !!----------------------------------------------------------------------
-      USE oce, ONLY:   tsa                         ! (tsa) used as 2 3D workspace
-      !!
       INTEGER, INTENT(in) ::   kt    ! ocean time-step index
       !!
@@ -199,7 +194,6 @@
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj 
       !!----------------------------------------------------------------------
-      
-      zhpi => tsa(:,:,:,1) 
-      zhpj => tsa(:,:,:,2) 
+      !  
+      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj )
       !
       IF( kt == nit000 ) THEN
@@ -245,4 +239,6 @@
       END DO
       !
+      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj )
+      !
    END SUBROUTINE hpg_zco
 
@@ -256,6 +252,4 @@
       !! ** Action  : - Update (ua,va) with the now hydrastatic pressure trend
       !!---------------------------------------------------------------------- 
-      USE oce, ONLY:   tsa                         ! (tsa) used as 2 3D workspace
-      !!
       INTEGER, INTENT(in) ::   kt    ! ocean time-step index
       !!
@@ -265,7 +259,6 @@
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj 
       !!----------------------------------------------------------------------
-       
-      zhpi => tsa(:,:,:,1) 
-      zhpj => tsa(:,:,:,2) 
+      !
+      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj )
       !
       IF( kt == nit000 ) THEN
@@ -343,5 +336,6 @@
       END DO
       !
-
+      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj )
+      !
    END SUBROUTINE hpg_zps
 
@@ -365,6 +359,4 @@
       !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
       !!----------------------------------------------------------------------
-      USE oce, ONLY:   tsa                         ! (tsa) used as 2 3D workspace
-      !!
       INTEGER, INTENT(in) ::   kt    ! ocean time-step index
       !!
@@ -373,7 +365,6 @@
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj 
       !!----------------------------------------------------------------------
-
-      zhpi => tsa(:,:,:,1) 
-      zhpj => tsa(:,:,:,2) 
+      !
+      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj )
       !
       IF( kt == nit000 ) THEN
@@ -432,4 +423,6 @@
       END DO
       !
+      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj )
+      !
    END SUBROUTINE hpg_sco
 
@@ -442,6 +435,4 @@
       !! Reference: Shchepetkin and McWilliams, J. Geophys. Res., 108(C3), 3090, 2003
       !!----------------------------------------------------------------------
-      USE oce     , ONLY:   tsa                         ! (tsa) used as 2 3D workspace
-      !!
       INTEGER, INTENT(in) ::   kt    ! ocean time-step index
       !!
@@ -458,7 +449,5 @@
       CALL wrk_alloc( jpi, jpj, jpk, dzx  , dzy  , dzz  , dzu  , dzv  , dzw   ) 
       CALL wrk_alloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 
-      CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k                      ) 
-      zhpi => tsa(:,:,:,1) 
-      zhpj => tsa(:,:,:,2) 
+      CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k,  zhpi,  zhpj        ) 
       !
 
@@ -660,5 +649,5 @@
       CALL wrk_dealloc( jpi, jpj, jpk, dzx  , dzy  , dzz  , dzu  , dzv  , dzw   ) 
       CALL wrk_dealloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 
-      CALL wrk_dealloc( jpi, jpj, jpk, rho_i, rho_j, rho_k                      ) 
+      CALL wrk_dealloc( jpi, jpj, jpk, rho_i, rho_j, rho_k,  zhpi,  zhpj        ) 
       !
    END SUBROUTINE hpg_djc
@@ -678,7 +667,4 @@
       !!
       !!----------------------------------------------------------------------
-      USE oce     , ONLY:   tsa                          ! (tsa) used as 2 3D workspace
-      !!----------------------------------------------------------------------
-      !!
       INTEGER, PARAMETER  :: polynomial_type = 1    ! 1: cubic spline, 2: linear
       INTEGER, INTENT(in) ::   kt                   ! ocean time-step index
@@ -697,7 +683,6 @@
       !!----------------------------------------------------------------------
       !
-      CALL wrk_alloc( jpi, jpj, jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 
-      zdeptht => tsa(:,:,:,1) 
-      zrhh    => tsa(:,:,:,2) 
+      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 
+      CALL wrk_alloc( jpi,jpj,jpk, zdeptht, zrhh ) 
       !
       IF( kt == nit000 ) THEN
@@ -717,5 +702,5 @@
       zrhh(:,:,:) = rhd(:,:,:)
       
-      ! Preparing vertical density profile for hybrid-sco coordinate
+      ! Preparing vertical density profile "zrhh(:,:,:)" for hybrid-sco coordinate
       DO jj = 1, jpj
         DO ji = 1, jpi   
@@ -732,4 +717,5 @@
       END DO
 
+      ! Transfer the depth of "T(:,:,:)" to vertical coordinate "zdeptht(:,:,:)"
       DO jj = 1, jpj
         DO ji = 1, jpi
@@ -753,7 +739,8 @@
       ! Construct the vertical density profile with the 
       ! constrained cubic spline interpolation
+      ! rho(z) = asp + bsp*z + csp*z^2 + dsp*z^3
       CALL cspline(fsp,xsp,asp,bsp,csp,dsp,polynomial_type)      
 
-      ! Calculate the hydrostatic pressure at T(ji,jj,1)
+      ! Integrate the hydrostatic pressure "zhpi(:,:,:)" at "T(ji,jj,1)"
       DO jj = 2, jpj
         DO ji = 2, jpi 
@@ -768,5 +755,5 @@
       END DO
 
-      ! Calculate the pressure at T(ji,jj,2:jpkm1)
+      ! Calculate the pressure "zhpi(:,:,:)" at "T(ji,jj,2:jpkm1)"
       DO jk = 2, jpkm1                                  
         DO jj = 2, jpj     
@@ -953,5 +940,6 @@
       END DO
       !
-      CALL wrk_dealloc( jpi, jpj, jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 
+      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 
+      CALL wrk_dealloc( jpi,jpj,jpk, zdeptht, zrhh ) 
       !
    END SUBROUTINE hpg_prj
@@ -964,6 +952,5 @@
       !!          
       !! ** Method  :   f(x) = asp + bsp*x + csp*x^2 + dsp*x^3 
-      !! Reference: K.W. Brodlie, A review of mehtods for curve and function
-      !!                          drawing, 1980
+      !! Reference: CJC Kruger, Constrained Cubic Spline Interpoltation
       !!
       !!----------------------------------------------------------------------
@@ -987,72 +974,82 @@
       jpkm1 = size(fsp,3) - 1
 
-      ! Clean output arrays
-      asp = 0.0_wp
-      bsp = 0.0_wp
-      csp = 0.0_wp
-      dsp = 0.0_wp
       
-      DO ji = 1, jpi
-        DO jj = 1, jpj
-          IF (polynomial_type == 1) THEN     ! Constrained Cubic Spline
-             DO jk = 2, jpkm1-1
-                zdxtmp1 = xsp(ji,jj,jk)   - xsp(ji,jj,jk-1)  
-                zdxtmp2 = xsp(ji,jj,jk+1) - xsp(ji,jj,jk)  
-                zdf1    = ( fsp(ji,jj,jk)   - fsp(ji,jj,jk-1) ) / zdxtmp1
-                zdf2    = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk)   ) / zdxtmp2
+      IF (polynomial_type == 1) THEN     ! Constrained Cubic Spline
+         DO ji = 1, jpi
+            DO jj = 1, jpj
+           !!Fritsch&Butland's method, 1984 (preferred, but more computation)              
+           !    DO jk = 2, jpkm1-1
+           !       zdxtmp1 = xsp(ji,jj,jk)   - xsp(ji,jj,jk-1)  
+           !       zdxtmp2 = xsp(ji,jj,jk+1) - xsp(ji,jj,jk)  
+           !       zdf1    = ( fsp(ji,jj,jk)   - fsp(ji,jj,jk-1) ) / zdxtmp1
+           !       zdf2    = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk)   ) / zdxtmp2
+           !
+           !       zalpha = ( zdxtmp1 + 2._wp * zdxtmp2 ) / ( zdxtmp1 + zdxtmp2 ) / 3._wp
+           !     
+           !       IF(zdf1 * zdf2 <= 0._wp) THEN
+           !           zdf(jk) = 0._wp
+           !       ELSE
+           !         zdf(jk) = zdf1 * zdf2 / ( ( 1._wp - zalpha ) * zdf1 + zalpha * zdf2 )
+           !       ENDIF
+           !    END DO
+           
+           !!Simply geometric average
+               DO jk = 2, jpkm1-1
+                  zdf1 = (fsp(ji,jj,jk) - fsp(ji,jj,jk-1)) / (xsp(ji,jj,jk) - xsp(ji,jj,jk-1))
+                  zdf2 = (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / (xsp(ji,jj,jk+1) - xsp(ji,jj,jk))
+            
+                  IF(zdf1 * zdf2 <= 0._wp) THEN
+                     zdf(jk) = 0._wp
+                  ELSE
+                     zdf(jk) = 2._wp * zdf1 * zdf2 / (zdf1 + zdf2)
+                  ENDIF
+               END DO
+           
+               zdf(1)     = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / &
+                          &          ( xsp(ji,jj,2) - xsp(ji,jj,1) ) -  0.5_wp * zdf(2)
+               zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / &
+                          &          ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - &
+                          & 0.5_wp * zdf(jpkm1 - 1)
    
-                zalpha = ( zdxtmp1 + 2._wp * zdxtmp2 ) / ( zdxtmp1 + zdxtmp2 ) / 3._wp
-              
-                IF(zdf1 * zdf2 <= 0._wp) THEN
-                    zdf(jk) = 0._wp
-                ELSE
-                  zdf(jk) = zdf1 * zdf2 / ( ( 1._wp - zalpha ) * zdf1 + zalpha * zdf2 )
-                ENDIF
-             END DO
-
-             zdf(1)     = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / &
-                        &          ( xsp(ji,jj,2) - xsp(ji,jj,1) ) -  0.5_wp * zdf(2)
-             zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / &
-                        &          ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - &
-                        & 0.5_wp * zdf(jpkm1 - 1)
+               DO jk = 1, jpkm1 - 1
+                 zdxtmp = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 
+                 ztmp1  = (zdf(jk+1) + 2._wp * zdf(jk)) / zdxtmp
+                 ztmp2  =  6._wp * (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / zdxtmp / zdxtmp
+                 zddf1  = -2._wp * ztmp1 + ztmp2 
+                 ztmp1  = (2._wp * zdf(jk+1) + zdf(jk)) / zdxtmp
+                 zddf2  =  2._wp * ztmp1 - ztmp2 
+      
+                 dsp(ji,jj,jk) = (zddf2 - zddf1) / 6._wp / zdxtmp
+                 csp(ji,jj,jk) = ( xsp(ji,jj,jk+1) * zddf1 - xsp(ji,jj,jk)*zddf2 ) / 2._wp / zdxtmp
+                 bsp(ji,jj,jk) = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp - & 
+                               & csp(ji,jj,jk) * ( xsp(ji,jj,jk+1) + xsp(ji,jj,jk) ) - &
+                               & dsp(ji,jj,jk) * ((xsp(ji,jj,jk+1) + xsp(ji,jj,jk))**2 - &
+                               &                   xsp(ji,jj,jk+1) * xsp(ji,jj,jk))
+                 asp(ji,jj,jk) = fsp(ji,jj,jk) - xsp(ji,jj,jk) * (bsp(ji,jj,jk) + &
+                               &                (xsp(ji,jj,jk) * (csp(ji,jj,jk) + &
+                               &                 dsp(ji,jj,jk) * xsp(ji,jj,jk))))
+               END DO
+            END DO
+         END DO
+ 
+      ELSE IF (polynomial_type == 2) THEN     ! Linear
+         DO ji = 1, jpi
+            DO jj = 1, jpj
+               DO jk = 1, jpkm1-1
+                  zdxtmp =xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 
+                  ztmp1 = fsp(ji,jj,jk+1) - fsp(ji,jj,jk)
    
-             DO jk = 1, jpkm1 - 1
-                zdxtmp = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 
-                ztmp1  = (zdf(jk+1) + 2._wp * zdf(jk)) / zdxtmp
-                ztmp2  =  6._wp * (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / zdxtmp / zdxtmp
-                zddf1  = -2._wp * ztmp1 + ztmp2 
-                ztmp1  = (2._wp * zdf(jk+1) + zdf(jk)) / zdxtmp
-                zddf2  =  2._wp * ztmp1 - ztmp2 
-   
-                dsp(ji,jj,jk) = (zddf2 - zddf1) / 6._wp / zdxtmp
-                csp(ji,jj,jk) = ( xsp(ji,jj,jk+1) * zddf1 - xsp(ji,jj,jk)*zddf2 ) / 2._wp / zdxtmp
-                bsp(ji,jj,jk) = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp - & 
-                              & csp(ji,jj,jk) * ( xsp(ji,jj,jk+1) + xsp(ji,jj,jk) ) - &
-                              & dsp(ji,jj,jk) * ( xsp(ji,jj,jk+1)**2 + &
-                              &                   xsp(ji,jj,jk+1) * xsp(ji,jj,jk) + &
-                              &                   xsp(ji,jj,jk)**2 )
-                asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk) - &
-                              &                 csp(ji,jj,jk) * xsp(ji,jj,jk)**2 - &
-                              &                 dsp(ji,jj,jk) * xsp(ji,jj,jk)**3
-             END DO
- 
-          ELSE IF (polynomial_type == 2) THEN     ! Linear
- 
-             DO jk = 1, jpkm1-1
-                zdxtmp =xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 
-                ztmp1 = fsp(ji,jj,jk+1) - fsp(ji,jj,jk)
-   
-                dsp(ji,jj,jk) = 0._wp
-                csp(ji,jj,jk) = 0._wp
-                bsp(ji,jj,jk) = ztmp1 / zdxtmp
-                asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk)
-             END DO
-
-          ELSE
-             CALL ctl_stop( 'invalid polynomial type in cspline' )
-          ENDIF
-
-        END DO
-      END DO
+                  dsp(ji,jj,jk) = 0._wp
+                  csp(ji,jj,jk) = 0._wp
+                  bsp(ji,jj,jk) = ztmp1 / zdxtmp
+                  asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+
+      ELSE
+           CALL ctl_stop( 'invalid polynomial type in cspline' )
+      ENDIF
+
       
    END SUBROUTINE cspline
@@ -1069,5 +1066,4 @@
       !!                extrapolation is also permitted (no value limit) 
       !!
-      !! H.Liu, Jan 2009,  POL 
       !!----------------------------------------------------------------------
       IMPLICIT NONE
