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 3080)
+++ branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 3081)
@@ -94,5 +94,5 @@
       ENDIF      
       !
-      SELECT CASE ( nhpg )      ! Hydrastatic pressure gradient computation
+      SELECT CASE ( nhpg )      ! Hydrostatic pressure gradient computation
       CASE (  0 )   ;   CALL hpg_zco    ( kt )      ! z-coordinate
       CASE (  1 )   ;   CALL hpg_zps    ( kt )      ! z-coordinate plus partial steps (interpolation)
@@ -171,5 +171,5 @@
       IF( ln_hpg_prj )   nhpg = 7
       !
-      !                               ! Consitency check
+      !                               ! Consistency check
       ioptio = 0 
       IF( ln_hpg_zco )   ioptio = ioptio + 1
@@ -1045,14 +1045,13 @@
       !!
       INTEGER, PARAMETER  :: polynomial_type = 1    ! 1: cubic spline, 2: linear
-      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
-      !!
-      INTEGER  ::   ji, jj, jk, jkk                ! dummy loop indices
-      REAL(wp) ::   zcoef0, znad               ! temporary scalars
-      !!
-      !! The local varialbes for the correction term
-      INTEGER  :: jke, jkw, jkn, jks, jk1, jis, jid, jjs, jjd
-      REAL(wp) :: zze, zzw, zzn, zzs, rre, rrw, rrn, rrs
-      REAL(wp) :: zuijk, zvijk, pe, pw, pwes, pwed, pn, ps, pnss, pnsd, deps
-      REAL(wp) :: rhde, rhdw, rhdn, rhds,rhdt1, rhdt2
+      INTEGER, INTENT(in) ::   kt                   ! ocean time-step index
+      !!
+      INTEGER  ::   ji, jj, jk, jkk                 ! dummy loop indices
+      REAL(wp) ::   zcoef0, znad                    ! temporary scalars
+      !!
+      !! 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
@@ -1073,5 +1072,5 @@
       zcoef0 = - grav 
       znad = 0.0_wp
-      IF(lk_vvl) znad = 1._wp
+      IF( lk_vvl ) znad = 1._wp
 
       ! Save fsde3w and rhd
@@ -1079,20 +1078,19 @@
       rhd_tmp(:,:,:)    = rhd(:,:,:)   
 
-      ! Clean 3-D work arraies
+      ! Clean 3-D work arrays
       zhpi(:,:,:) = 0.
       
-      !how to add vector opt.? N.B., jpi&jpi rather than jpim1&jpjm1 are needed here
-
+      ! 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)
-            ELSE IF(jk < jpkm1) THEN
-            DO jkk = jk+1, jpk
-              rhd(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 
+          IF( jk <= 0 ) THEN; rhd(ji,jj,:) = 0._wp
+          ELSE IF(jk == 1) THEN; rhd(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), &
+                                         fsde3w(ji,jj,jkk-2), rhd(ji,jj,jkk-1), rhd(ji,jj,jkk-2))
+             END DO 
           END IF
         END DO
@@ -1104,5 +1102,5 @@
           fsde3w(ji,jj,1) = fsde3w(ji,jj,1) - sshn(ji,jj) * znad
           DO jk = 2, jpk
-            fsde3w(ji,jj,jk) = fsde3w(ji,jj,jk-1) + fse3w(ji,jj,jk)
+             fsde3w(ji,jj,jk) = fsde3w(ji,jj,jk-1) + fse3w(ji,jj,jk)
           END DO
         END DO
@@ -1118,6 +1116,6 @@
       END DO
 
-      !                 ! Construct the vertical density profile with the 
-      !                 !Constrained cubic spline interpolation
+      ! Construct the vertical density profile with the 
+      ! constrained cubic spline interpolation
       CALL cspline(fsp,xsp,asp,bsp,csp,dsp,polynomial_type)      
 
@@ -1126,9 +1124,9 @@
         DO ji = 2, jpi 
           rhdt1 = rhd(ji,jj,1) - interp3(fsde3w(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)
+                                         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
 
-     !                  ! assuming linear profile across the top half surface layer
+          ! assuming linear profile across the top half surface layer
           zhpi(ji,jj,1) =  0.5_wp * fse3w(ji,jj,1) * rhdt1  
         END DO
@@ -1139,8 +1137,8 @@
         DO jj = 2, jpj     
           DO ji = 2, jpi
-            zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + &
-                             integ2(fsde3w(ji,jj,jk-1),fsde3w(ji,jj,jk),&
-                                    asp(ji,jj,jk-1),bsp(ji,jj,jk-1),&
-                                    csp(ji,jj,jk-1),dsp(ji,jj,jk-1))
+            zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) +                          &
+                             integ2(fsde3w(ji,jj,jk-1), fsde3w(ji,jj,jk),&
+                                    asp(ji,jj,jk-1),    bsp(ji,jj,jk-1), &
+                                    csp(ji,jj,jk-1),    dsp(ji,jj,jk-1))
           END DO
         END DO
@@ -1148,9 +1146,8 @@
 
       ! Z coordinate of U(ji,jj,1:jpkm1) and V(ji,jj,1:jpkm1)
-
       DO jj = 2, jpjm1     
         DO ji = 2, jpim1  
-          zu(ji,jj,1) = - ( 0.5_wp * fse3uw(ji,jj,1) - sshu_n(ji,jj) * znad)
-          zv(ji,jj,1) = - ( 0.5_wp * fse3vw(ji,jj,1) - sshv_n(ji,jj) * znad)
+          zu(ji,jj,1) = - ( fse3u(ji,jj,1) - sshu_n(ji,jj) * znad)
+          zv(ji,jj,1) = - ( fse3v(ji,jj,1) - sshv_n(ji,jj) * znad)
         END DO
       END DO
@@ -1159,11 +1156,18 @@
         DO jj = 2, jpjm1     
           DO ji = 2, jpim1  
-            zu(ji,jj,jk) = zu(ji,jj,jk-1)- fse3uw(ji,jj,jk)
-            zv(ji,jj,jk) = zv(ji,jj,jk-1)- fse3vw(ji,jj,jk)
+            zu(ji,jj,jk) = zu(ji,jj,jk-1)- fse3u(ji,jj,jk)
+            zv(ji,jj,jk) = zv(ji,jj,jk-1)- fse3v(ji,jj,jk)
           END DO
         END DO
       END DO
                
-      !  Start pressure integration
+      DO jk = 1, jpkm1                                  
+        DO jj = 2, jpjm1     
+          DO ji = 2, jpim1  
+            zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * fse3u(ji,jj,jk)
+            zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * fse3v(ji,jj,jk)
+          END DO
+        END DO
+      END DO
 
       DO jk = 1, jpkm1                                  
@@ -1176,171 +1180,134 @@
 
             !!!!!     for u equation
-            IF( -fsde3w(ji+1,jj,mbathy(ji+1,jj)) >= -fsde3w(ji,jj,mbathy(ji,jj)) ) THEN
-              jis = ji + 1; jid = ji
-            ELSE
-              jis = ji;     jid = ji +1
+            IF( jk <= mbku(ji,jj) ) THEN
+               IF( -fsde3w(ji+1,jj,mbku(ji,jj)) >= -fsde3w(ji,jj,mbku(ji,jj)) ) THEN
+                 jis = ji + 1; jid = ji
+               ELSE
+                 jis = ji;     jid = ji +1
+               ENDIF
+
+               ! integrate the pressure on the shallow side
+               jk1 = jk 
+               bhitwe = 0
+               DO WHILE ( -fsde3w(jis,jj,jk1) > zuijk )
+                 IF( jk1 == mbku(ji,jj) ) THEN
+                   bhitwe = 1
+                   EXIT
+                 ENDIF
+                 deps = min(fsde3w(jis,jj,jk1+1), -zuijk)
+                 pwes = pwes +                                    & 
+                      integ2(fsde3w(jis,jj,jk1), deps,            &
+                             asp(jis,jj,jk1),    bsp(jis,jj,jk1), &
+                             csp(jis,jj,jk1),    dsp(jis,jj,jk1))
+                 jk1 = jk1 + 1
+               END DO
+            
+               IF(bhitwe == 1) THEN
+                 zuijk = -fsde3w(jis,jj,jk1)
+               ENDIF
+
+               ! integrate the pressure on the deep side
+               jk1 = jk 
+               bhitwe = 0
+               DO WHILE ( -fsde3w(jid,jj,jk1) < zuijk )
+                 IF( jk1 == 1 ) THEN
+                   bhitwe = 1
+                   EXIT
+                 ENDIF
+                 deps = max(fsde3w(jid,jj,jk1-1), -zuijk)
+                 pwed = pwed +                                        & 
+                        integ2(deps,              fsde3w(jid,jj,jk1), &
+                               asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1),  &
+                               csp(jid,jj,jk1-1), dsp(jid,jj,jk1-1) )
+                 jk1 = jk1 - 1
+               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), &
+                                                 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
+               ENDIF
+
+
+               dpdx1 = 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)) ) 
+                ELSE
+                 dpdx2 = zcoef0 / e1u(ji,jj) * REAL(jis-jid, wp) * (pwes + pwed) 
+               ENDIF
+
+               ua(ji,jj,jk) = ua(ji,jj,jk) + (dpdx1 + dpdx2) * &
+               &           umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk)
             ENDIF
-
-            !     !integrate the pressure on the shallow side
-            jk1 = jk 
-            bhitwe = 0
-            IF( jk1 == mbathy(jis,jj) ) THEN
-              bhitwe = 1
-            ELSE
-              DO WHILE ( -fsde3w(jis,jj,jk1+1) > zuijk )
-                pwes = pwes + & 
-                       integ2(fsde3w(jis,jj,jk1), fsde3w(jis,jj,jk1+1),&
-                              asp(jis,jj,jk1)   , bsp(jis,jj,jk1),     &
-                              csp(jis,jj,jk1)   , dsp(jis,jj,jk1))
-                jk1 = jk1 + 1
-                IF( jk1 == mbathy(jis,jj) ) THEN
-                  bhitwe = 1
-                  EXIT
-                END IF
-              END DO
-            ENDIF
-            
-            IF( bhitwe /= 1 ) THEN
-              pwes = pwes + & 
-                     integ2(fsde3w(jis,jj,jk1), -zuijk,      &
-                            asp(jis,jj,jk1), bsp(jis,jj,jk1),&
-                            csp(jis,jj,jk1), dsp(jis,jj,jk1))
-            ELSE
-               zuijk = -fsde3w(jis,jj,jk1)
-            ENDIF
-
-            !     !integrate the pressure on the deep side
-            jk1 = jk 
-            bhitwe = 0
-            IF( jk1 == 1 ) THEN
-              bhitwe = 1
-            ELSE
-              DO WHILE ( -fsde3w(jid,jj,jk1-1) < zuijk )
-                pwed = pwed + & 
-                       integ2(fsde3w(jid,jj,jk1-1), fsde3w(jid,jj,jk1),&
-                              asp(jid,jj,jk1-1),    bsp(jid,jj,jk1-1), &
-                              csp(jid,jj,jk1-1),    dsp(jid,jj,jk1-1))
-                jk1 = jk1 - 1
-                IF( jk1 == 1 ) THEN
-                  bhitwe = 1
-                  EXIT
-                END IF
-              END DO
-            ENDIF
-            
-            IF( bhitwe /= 1 ) THEN
-              pwed = pwed + & 
-                     integ2(-zuijk,             fsde3w(jid,jj,jk1),&
-                             asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1), &
-                             csp(jid,jj,jk1-1), dsp(jid,jj,jk1-1))
-            ELSE
-              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), &
-                                              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
-            ENDIF
-
-            IF( jid > jis ) THEN
-              pe = pwed; pw = pwes
-            ELSE
-              pe = pwes; pw = pwed
-            ENDIF
-
-            dpdx1 = zcoef0 / e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk))
-            IF( lk_vvl ) THEN
-              dpdx2 = zcoef0 / e1u(ji,jj) * (pw + pe + (sshn(ji+1,jj)-sshn(ji,jj))) 
-             ELSE
-              dpdx2 = zcoef0 / e1u(ji,jj) * (pw + pe) 
-            ENDIF
-
-            ua(ji,jj,jk) = ua(ji,jj,jk) + (dpdx1 + dpdx2) * &
-               &           umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk)
   
             !!!!!     for v equation
-
-            IF( -fsde3w(ji,jj+1,mbathy(ji,jj+1)) >= -fsde3w(ji,jj,mbathy(ji,jj)) ) THEN
-              jjs = jj + 1; jjd = jj
-            ELSE
-              jjs = jj    ; jjd = jj + 1
+            IF( jk <= mbkv(ji,jj) ) THEN
+               IF( -fsde3w(ji,jj+1,mbkv(ji,jj)) >= -fsde3w(ji,jj,mbkv(ji,jj)) ) THEN
+                 jjs = jj + 1; jjd = jj
+               ELSE
+                 jjs = jj    ; jjd = jj + 1
+               ENDIF
+
+               ! integrate the pressure on the shallow side
+               jk1 = jk 
+               bhitns = 0
+               DO WHILE ( -fsde3w(ji,jjs,jk1) > zvijk )
+                 IF( jk1 == mbkv(ji,jj) ) THEN
+                   bhitns = 1
+                   EXIT
+                 ENDIF
+                 deps = min(fsde3w(ji,jjs,jk1+1), -zvijk)
+                 pnss = pnss +                                      & 
+                        integ2(fsde3w(ji,jjs,jk1), deps,            &
+                               asp(ji,jjs,jk1),    bsp(ji,jjs,jk1), &
+                               csp(ji,jjs,jk1),    dsp(ji,jjs,jk1) )
+                 jk1 = jk1 + 1
+               END DO
+            
+               IF(bhitns == 1) THEN
+                 zvijk = -fsde3w(ji,jjs,jk1)
+               ENDIF
+
+               ! integrate the pressure on the deep side
+               jk1 = jk 
+               bhitns = 0
+               DO WHILE ( -fsde3w(ji,jjd,jk1) < zvijk )
+                 IF( jk1 == 1 ) THEN
+                   bhitns = 1
+                   EXIT
+                 ENDIF
+                 deps = max(fsde3w(ji,jjd,jk1-1), -zvijk)
+                 pnsd = pnsd +                                        & 
+                        integ2(deps,              fsde3w(ji,jjd,jk1), &
+                               asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1), &
+                               csp(ji,jjd,jk1-1), dsp(ji,jjd,jk1-1) )
+                 jk1 = jk1 - 1
+               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), &
+                                                 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
+               ENDIF
+
+
+               dpdy1 = 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)) ) 
+               ELSE
+                   dpdy2 = zcoef0 / e2v(ji,jj) * REAL(jjs-jjd, wp) * (pnss + pnsd ) 
+               ENDIF
+
+               va(ji,jj,jk) = va(ji,jj,jk) + (dpdy1 + dpdy2)*&
+               &              vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk)
             ENDIF
-
-            !     !integrate the pressure on the shallow side
-            jk1 = jk 
-            bhitns = 0
-            IF( jk1 == mbathy(ji,jjs) ) THEN
-              bhitns = 1
-            ELSE
-              DO WHILE ( -fsde3w(ji,jjs,jk1+1) > zvijk )
-                pnss = pnss + & 
-                       integ2(fsde3w(ji,jjs,jk1), fsde3w(ji,jjs,jk1+1),&
-                              asp(ji,jjs,jk1),    bsp(ji,jjs,jk1),     &
-                              csp(ji,jjs,jk1),    dsp(ji,jjs,jk1) )
-                jk1 = jk1 + 1
-                IF( jk1 == mbathy(ji,jjs) ) THEN
-                  bhitns = 1
-                  EXIT
-                END IF
-              END DO
-            ENDIF
-            
-            IF( bhitns /= 1 ) THEN
-              pnss = pnss + & 
-                     integ2(fsde3w(ji,jjs,jk1), -zvijk,          &
-                            asp(ji,jjs,jk1),    bsp(ji,jjs,jk1), &
-                            csp(ji,jjs,jk1),    dsp(ji,jjs,jk1) )
-            ELSE
-               zvijk = -fsde3w(ji,jjs,jk1)
-            ENDIF
-
-            !     !integrate the pressure on the deep side
-            jk1 = jk 
-            bhitns = 0
-            IF( jk1 == 1 ) THEN
-              bhitns = 1
-            ELSE
-              DO WHILE ( -fsde3w(ji,jjd,jk1-1) < zvijk )
-                pnsd = pnsd + & 
-                       integ2(fsde3w(ji,jjd,jk1-1), fsde3w(ji,jjd,jk1), &
-                              asp(ji,jjd,jk1-1),    bsp(ji,jjd,jk1-1),  &
-                              csp(ji,jjd,jk1-1),    dsp(ji,jjd,jk1-1) )
-                jk1 = jk1 - 1
-                IF( jk1 == 1 ) THEN
-                  bhitns = 1
-                  EXIT
-                END IF
-              END DO
-            ENDIF
-            
-            IF( bhitns /= 1 ) THEN
-              pnsd = pnsd + & 
-                     integ2(-zvijk,            fsde3w(ji,jjd,jk1), &
-                            asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1),  &
-                            csp(ji,jjd,jk1-1), dsp(ji,jjd,jk1-1) )
-            ELSE
-              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), &
-                                              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
-            ENDIF
-
-            IF( jjd > jjs ) THEN
-              pn = pnsd; ps = pnss
-            ELSE
-              pn = pnss; ps = pnsd
-            ENDIF
-
-            dpdy1 = zcoef0 / e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk))
-            if( lk_vvl ) then
-                dpdy2 = zcoef0 / e2v(ji,jj) * (ps + pn + (sshn(ji,jj+1)-sshn(ji,jj))) 
-              else
-                dpdy2 = zcoef0 / e2v(ji,jj) * (ps + pn ) 
-            end if
-
-            va(ji,jj,jk) = va(ji,jj,jk) + (dpdy1 + dpdy2)* &
-            &              vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk)
 
                     
@@ -1399,5 +1366,5 @@
       Do jj = 1, jpj
 
-      If (polynomial_type == 1) Then     !Constrained Cubic Spline
+      If (polynomial_type == 1) Then     ! Constrained Cubic Spline
           Do jk = 2, jpkm1-1
              zdxtmp1 = xsp(ji,jj,jk)   - xsp(ji,jj,jk-1)  
@@ -1441,5 +1408,5 @@
           End Do
 
-       Else If (polynomial_type == 2) Then     !Linear
+       Else If (polynomial_type == 2) Then     ! Linear
 
         Do jk = 1, jpkm1-1
