Index: branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 3081)
+++ branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 3099)
@@ -1020,7 +1020,7 @@
       !!
       !! ** Method  :   s-coordinate case.
-      !!      Reformulate the horizontal hydrostatical pressure gradient
-      !!      term using Pressure Jacobian. A new correction term 
-      !!      is developed to eliminate the sigma-coordinate error.
+      !!      A Pressure-Jacobian horizontal pressure gradient method
+      !!      based on the constrained cubic-spline interpolation for
+      !!      all vertical coordinate systems
       !!
       !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
@@ -1030,6 +1030,6 @@
 
       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
-      USE oce     , ONLY:   fsde3w_tmp => ta       ! (ta,sa) used as 3D workspace
-      USE oce     , ONLY:      rhd_tmp => sa 
+      USE oce     , ONLY:       zdeptht => ta       ! (ta,sa) used as 3D workspace
+      USE oce     , ONLY:          zrhh => sa 
       USE wrk_nemo, ONLY:         zhpi => wrk_3d_3 
       USE wrk_nemo, ONLY:           zu => wrk_3d_4 
@@ -1052,8 +1052,8 @@
       !! The local variables for the correction term
       INTEGER  :: jk1, jis, jid, jjs, jjd
-      REAL(wp) :: zuijk, zvijk, pwes, pwed, pnss, pnsd, deps
-      REAL(wp) :: rhdt1 
-      REAL(wp) :: dpdx1, dpdx2, dpdy1, dpdy2
-      INTEGER  :: bhitwe, bhitns
+      REAL(wp) :: zuijk, zvijk, zpwes, zpwed, zpnss, zpnsd, zdeps
+      REAL(wp) :: zrhdt1 
+      REAL(wp) :: zdpdx1, zdpdx2, zdpdy1, zdpdy2
+      INTEGER  :: zbhitwe, zbhitns
       !!----------------------------------------------------------------------
 
@@ -1074,24 +1074,19 @@
       IF( lk_vvl ) znad = 1._wp
 
-      ! Save fsde3w and rhd
-      fsde3w_tmp(:,:,:) = fsde3w(:,:,:) 
-      rhd_tmp(:,:,:)    = rhd(:,:,:)   
-
       ! Clean 3-D work arrays
       zhpi(:,:,:) = 0.
       
-      ! how to add vector opt.? N.B., jpi&jpi rather than jpim1&jpjm1 are needed here
       ! Preparing vertical density profile for hybrid-sco coordinate
       DO jj = 1, jpj
         DO ji = 1, jpi   
           jk = mbathy(ji,jj)
-          IF( jk <= 0 ) THEN; rhd(ji,jj,:) = 0._wp
-          ELSE IF(jk == 1) THEN; rhd(ji,jj, jk+1:jpk) = rhd(ji,jj,jk)
+          IF( jk <= 0 ) THEN; zrhh(ji,jj,:) = 0._wp
+          ELSE IF(jk == 1) THEN; zrhh(ji,jj, jk+1:jpk) = rhd(ji,jj,jk)
           ELSE IF(jk < jpkm1) THEN
              DO jkk = jk+1, jpk
-                rhd(ji,jj,jkk) = interp1(fsde3w(ji,jj,jkk),   fsde3w(ji,jj,jkk-1), &
+                zrhh(ji,jj,jkk) = interp1(fsde3w(ji,jj,jkk),   fsde3w(ji,jj,jkk-1), &
                                          fsde3w(ji,jj,jkk-2), rhd(ji,jj,jkk-1), rhd(ji,jj,jkk-2))
              END DO 
-          END IF
+          ENDIF
         END DO
       END DO
@@ -1099,8 +1094,8 @@
       DO jj = 1, jpj
         DO ji = 1, jpi
-          fsde3w(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1)
-          fsde3w(ji,jj,1) = fsde3w(ji,jj,1) - sshn(ji,jj) * znad
+          zdeptht(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1)
+          zdeptht(ji,jj,1) = zdeptht(ji,jj,1) - sshn(ji,jj) * znad
           DO jk = 2, jpk
-             fsde3w(ji,jj,jk) = fsde3w(ji,jj,jk-1) + fse3w(ji,jj,jk)
+             zdeptht(ji,jj,jk) = zdeptht(ji,jj,jk-1) + fse3w(ji,jj,jk)
           END DO
         END DO
@@ -1110,6 +1105,6 @@
         DO jj = 1, jpj
           DO ji = 1, jpi
-            fsp(ji,jj,jk) = rhd(ji,jj,jk)
-            xsp(ji,jj,jk) = fsde3w(ji,jj,jk)
+            fsp(ji,jj,jk) = zrhh(ji,jj,jk)
+            xsp(ji,jj,jk) = zdeptht(ji,jj,jk)
           END DO
         END DO
@@ -1123,11 +1118,11 @@
       DO jj = 2, jpj
         DO ji = 2, jpi 
-          rhdt1 = rhd(ji,jj,1) - interp3(fsde3w(ji,jj,1),asp(ji,jj,1), &
+          zrhdt1 = zrhh(ji,jj,1) - interp3(zdeptht(ji,jj,1),asp(ji,jj,1), &
                                          bsp(ji,jj,1),   csp(ji,jj,1), &
-                                         dsp(ji,jj,1) ) * 0.5_wp * fsde3w(ji,jj,1)
-          rhdt1 = max(rhdt1, 1000._wp - rau0)        ! no lighter than fresh water
+                                         dsp(ji,jj,1) ) * 0.5_wp * zdeptht(ji,jj,1)
+          zrhdt1 = MAX(zrhdt1, 1000._wp - rau0)        ! no lighter than fresh water
 
           ! assuming linear profile across the top half surface layer
-          zhpi(ji,jj,1) =  0.5_wp * fse3w(ji,jj,1) * rhdt1  
+          zhpi(ji,jj,1) =  0.5_wp * fse3w(ji,jj,1) * zrhdt1  
         END DO
       END DO
@@ -1138,5 +1133,5 @@
           DO ji = 2, jpi
             zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) +                          &
-                             integ2(fsde3w(ji,jj,jk-1), fsde3w(ji,jj,jk),&
+                             integ2(zdeptht(ji,jj,jk-1), zdeptht(ji,jj,jk),&
                                     asp(ji,jj,jk-1),    bsp(ji,jj,jk-1), &
                                     csp(ji,jj,jk-1),    dsp(ji,jj,jk-1))
@@ -1174,6 +1169,6 @@
         DO jj = 2, jpjm1     
           DO ji = 2, jpim1  
-            pwes = 0._wp; pwed = 0._wp
-            pnss = 0._wp; pnsd = 0._wp
+            zpwes = 0._wp; zpwed = 0._wp
+            zpnss = 0._wp; zpnsd = 0._wp
             zuijk = zu(ji,jj,jk)
             zvijk = zv(ji,jj,jk)
@@ -1181,5 +1176,5 @@
             !!!!!     for u equation
             IF( jk <= mbku(ji,jj) ) THEN
-               IF( -fsde3w(ji+1,jj,mbku(ji,jj)) >= -fsde3w(ji,jj,mbku(ji,jj)) ) THEN
+               IF( -zdeptht(ji+1,jj,mbku(ji,jj)) >= -zdeptht(ji,jj,mbku(ji,jj)) ) THEN
                  jis = ji + 1; jid = ji
                ELSE
@@ -1189,13 +1184,13 @@
                ! integrate the pressure on the shallow side
                jk1 = jk 
-               bhitwe = 0
-               DO WHILE ( -fsde3w(jis,jj,jk1) > zuijk )
+               zbhitwe = 0
+               DO WHILE ( -zdeptht(jis,jj,jk1) > zuijk )
                  IF( jk1 == mbku(ji,jj) ) THEN
-                   bhitwe = 1
+                   zbhitwe = 1
                    EXIT
                  ENDIF
-                 deps = min(fsde3w(jis,jj,jk1+1), -zuijk)
-                 pwes = pwes +                                    & 
-                      integ2(fsde3w(jis,jj,jk1), deps,            &
+                 zdeps = MIN(zdeptht(jis,jj,jk1+1), -zuijk)
+                 zpwes = zpwes +                                    & 
+                      integ2(zdeptht(jis,jj,jk1), zdeps,            &
                              asp(jis,jj,jk1),    bsp(jis,jj,jk1), &
                              csp(jis,jj,jk1),    dsp(jis,jj,jk1))
@@ -1203,19 +1198,19 @@
                END DO
             
-               IF(bhitwe == 1) THEN
-                 zuijk = -fsde3w(jis,jj,jk1)
+               IF(zbhitwe == 1) THEN
+                 zuijk = -zdeptht(jis,jj,jk1)
                ENDIF
 
                ! integrate the pressure on the deep side
                jk1 = jk 
-               bhitwe = 0
-               DO WHILE ( -fsde3w(jid,jj,jk1) < zuijk )
+               zbhitwe = 0
+               DO WHILE ( -zdeptht(jid,jj,jk1) < zuijk )
                  IF( jk1 == 1 ) THEN
-                   bhitwe = 1
+                   zbhitwe = 1
                    EXIT
                  ENDIF
-                 deps = max(fsde3w(jid,jj,jk1-1), -zuijk)
-                 pwed = pwed +                                        & 
-                        integ2(deps,              fsde3w(jid,jj,jk1), &
+                 zdeps = MAX(zdeptht(jid,jj,jk1-1), -zuijk)
+                 zpwed = zpwed +                                        & 
+                        integ2(zdeps,              zdeptht(jid,jj,jk1), &
                                asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1),  &
                                csp(jid,jj,jk1-1), dsp(jid,jj,jk1-1) )
@@ -1223,23 +1218,24 @@
                END DO
             
-               IF( bhitwe == 1 ) THEN
-                 deps = fsde3w(jid,jj,1) + min(zuijk, sshn(jid,jj)*znad)
-                 rhdt1 = rhd(jid,jj,1) - interp3(fsde3w(jid,jj,1), asp(jid,jj,1), &
+               IF( zbhitwe == 1 ) THEN
+                 zdeps = zdeptht(jid,jj,1) + MIN(zuijk, sshn(jid,jj)*znad)
+                 zrhdt1 = zrhh(jid,jj,1) - interp3(zdeptht(jid,jj,1), asp(jid,jj,1), &
                                                  bsp(jid,jj,1),    csp(jid,jj,1), &
-                                                 dsp(jid,jj,1)) * deps
-                 rhdt1 = max(rhdt1, 1000._wp - rau0)        ! no lighter than fresh water
-                 pwed  = pwed + 0.5_wp * (rhd(jid,jj,1) + rhdt1) * deps
+                                                 dsp(jid,jj,1)) * zdeps
+                 zrhdt1 = MAX(zrhdt1, 1000._wp - rau0)        ! no lighter than fresh water
+                 zpwed  = zpwed + 0.5_wp * (zrhh(jid,jj,1) + zrhdt1) * zdeps
                ENDIF
 
-
-               dpdx1 = zcoef0 / e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk))
+               ! update the momentum trends in u direction
+
+               zdpdx1 = zcoef0 / e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk))
                IF( lk_vvl ) THEN
-                 dpdx2 = zcoef0 / e1u(ji,jj) * & 
-                         ( REAL(jis-jid, wp) * (pwes + pwed) + (sshn(ji+1,jj)-sshn(ji,jj)) ) 
+                 zdpdx2 = zcoef0 / e1u(ji,jj) * & 
+                         ( REAL(jis-jid, wp) * (zpwes + zpwed) + (sshn(ji+1,jj)-sshn(ji,jj)) ) 
                 ELSE
-                 dpdx2 = zcoef0 / e1u(ji,jj) * REAL(jis-jid, wp) * (pwes + pwed) 
+                 zdpdx2 = zcoef0 / e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 
                ENDIF
 
-               ua(ji,jj,jk) = ua(ji,jj,jk) + (dpdx1 + dpdx2) * &
+               ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * &
                &           umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk)
             ENDIF
@@ -1247,5 +1243,5 @@
             !!!!!     for v equation
             IF( jk <= mbkv(ji,jj) ) THEN
-               IF( -fsde3w(ji,jj+1,mbkv(ji,jj)) >= -fsde3w(ji,jj,mbkv(ji,jj)) ) THEN
+               IF( -zdeptht(ji,jj+1,mbkv(ji,jj)) >= -zdeptht(ji,jj,mbkv(ji,jj)) ) THEN
                  jjs = jj + 1; jjd = jj
                ELSE
@@ -1255,13 +1251,13 @@
                ! integrate the pressure on the shallow side
                jk1 = jk 
-               bhitns = 0
-               DO WHILE ( -fsde3w(ji,jjs,jk1) > zvijk )
+               zbhitns = 0
+               DO WHILE ( -zdeptht(ji,jjs,jk1) > zvijk )
                  IF( jk1 == mbkv(ji,jj) ) THEN
-                   bhitns = 1
+                   zbhitns = 1
                    EXIT
                  ENDIF
-                 deps = min(fsde3w(ji,jjs,jk1+1), -zvijk)
-                 pnss = pnss +                                      & 
-                        integ2(fsde3w(ji,jjs,jk1), deps,            &
+                 zdeps = MIN(zdeptht(ji,jjs,jk1+1), -zvijk)
+                 zpnss = zpnss +                                      & 
+                        integ2(zdeptht(ji,jjs,jk1), zdeps,            &
                                asp(ji,jjs,jk1),    bsp(ji,jjs,jk1), &
                                csp(ji,jjs,jk1),    dsp(ji,jjs,jk1) )
@@ -1269,19 +1265,19 @@
                END DO
             
-               IF(bhitns == 1) THEN
-                 zvijk = -fsde3w(ji,jjs,jk1)
+               IF(zbhitns == 1) THEN
+                 zvijk = -zdeptht(ji,jjs,jk1)
                ENDIF
 
                ! integrate the pressure on the deep side
                jk1 = jk 
-               bhitns = 0
-               DO WHILE ( -fsde3w(ji,jjd,jk1) < zvijk )
+               zbhitns = 0
+               DO WHILE ( -zdeptht(ji,jjd,jk1) < zvijk )
                  IF( jk1 == 1 ) THEN
-                   bhitns = 1
+                   zbhitns = 1
                    EXIT
                  ENDIF
-                 deps = max(fsde3w(ji,jjd,jk1-1), -zvijk)
-                 pnsd = pnsd +                                        & 
-                        integ2(deps,              fsde3w(ji,jjd,jk1), &
+                 zdeps = MAX(zdeptht(ji,jjd,jk1-1), -zvijk)
+                 zpnsd = zpnsd +                                        & 
+                        integ2(zdeps,              zdeptht(ji,jjd,jk1), &
                                asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1), &
                                csp(ji,jjd,jk1-1), dsp(ji,jjd,jk1-1) )
@@ -1289,23 +1285,24 @@
                END DO
             
-               IF( bhitns == 1 ) THEN
-                 deps = fsde3w(ji,jjd,1) + min(zvijk, sshn(ji,jjd)*znad)
-                 rhdt1 = rhd(ji,jjd,1) - interp3(fsde3w(ji,jjd,1), asp(ji,jjd,1), &
+               IF( zbhitns == 1 ) THEN
+                 zdeps = zdeptht(ji,jjd,1) + MIN(zvijk, sshn(ji,jjd)*znad)
+                 zrhdt1 = zrhh(ji,jjd,1) - interp3(zdeptht(ji,jjd,1), asp(ji,jjd,1), &
                                                  bsp(ji,jjd,1),    csp(ji,jjd,1), &
-                                                 dsp(ji,jjd,1) ) * deps
-                 rhdt1 = max(rhdt1, 1000._wp - rau0)        ! no lighter than fresh water
-                 pnsd  = pnsd + 0.5_wp * (rhd(ji,jjd,1) + rhdt1) * deps
+                                                 dsp(ji,jjd,1) ) * zdeps
+                 zrhdt1 = MAX(zrhdt1, 1000._wp - rau0)        ! no lighter than fresh water
+                 zpnsd  = zpnsd + 0.5_wp * (zrhh(ji,jjd,1) + zrhdt1) * zdeps
                ENDIF
 
-
-               dpdy1 = zcoef0 / e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk))
+               ! update the momentum trends in v direction
+
+               zdpdy1 = zcoef0 / e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk))
                IF( lk_vvl ) THEN
-                   dpdy2 = zcoef0 / e2v(ji,jj) * &
-                           ( REAL(jjs-jjd, wp) * (pnss + pnsd) + (sshn(ji,jj+1)-sshn(ji,jj)) ) 
+                   zdpdy2 = zcoef0 / e2v(ji,jj) * &
+                           ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (sshn(ji,jj+1)-sshn(ji,jj)) ) 
                ELSE
-                   dpdy2 = zcoef0 / e2v(ji,jj) * REAL(jjs-jjd, wp) * (pnss + pnsd ) 
+                   zdpdy2 = zcoef0 / e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 
                ENDIF
 
-               va(ji,jj,jk) = va(ji,jj,jk) + (dpdy1 + dpdy2)*&
+               va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2)*&
                &              vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk)
             ENDIF
@@ -1316,12 +1313,9 @@
       END DO
 
-      ! Restore fsde3w and rhd
-      fsde3w(:,:,:) = fsde3w_tmp(:,:,:)
-      rhd(:,:,:)    = rhd_tmp(:,:,:)
-
       !
       IF( wrk_not_released(3, 3,4,5,6,7,8,9,10,11) )   &
          CALL ctl_stop('dyn:hpg_prj: failed to release workspace arrays')
       !
+XXXXXXX
    END SUBROUTINE hpg_prj
 
@@ -1362,70 +1356,68 @@
       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
+   
+                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)**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
       
-      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
-  
-             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 )
-             End If
-          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)**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' )
-      End If
-
-      End Do
-      End Do
-      
-   End Subroutine cspline
+   END SUBROUTINE cspline
 
 
@@ -1451,7 +1443,7 @@
       IF(abs(zdeltx) <= 10._wp * EPSILON(x)) THEN
         f = 0.5_wp * (fl + fr)
-       ELSE
+      ELSE
         f = ( (x - xl ) * fr - ( x - xr ) * fl ) / zdeltx
-      END IF
+      ENDIF
       
    END FUNCTION interp1
@@ -1480,7 +1472,8 @@
       !!                 ***  ROUTINE interp1  ***
       !!       
-      !! ** Purpose :   deriavtive of a cubic spline function
+      !! ** Purpose :   Calculate the first order of deriavtive of
+      !!                a cubic spline function y=a+b*x+c*x^2+d*x^3
       !!          
-      !! ** Method  :  
+      !! ** Method  :   f=dy/dx=b+2*c*x+3*d*x^2
       !!
       !!----------------------------------------------------------------------
