Index: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 5863)
+++ branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 5864)
@@ -234,8 +234,6 @@
             zcoef1 = zcoef0 * e3w_n(ji,jj,1)
             ! hydrostatic pressure gradient
-            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) / e1u(ji,jj)
-            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) / e2v(ji,jj)
-!!gm            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj)
-!!gm            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj)
+            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj)
+            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj)
             ! add to the general momentum trend
             ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
@@ -253,11 +251,9 @@
                zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   &
                   &           + zcoef1 * (  ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) )   &
-                  &                       - ( rhd(ji  ,jj,jk)+rhd(ji  ,jj,jk-1) )  ) / e1u(ji,jj)
-!!gm                  &                       - ( rhd(ji  ,jj,jk)+rhd(ji  ,jj,jk-1) )  ) * r1_e1u(ji,jj)
+                  &                       - ( rhd(ji  ,jj,jk)+rhd(ji  ,jj,jk-1) )  ) * r1_e1u(ji,jj)
 
                zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)   &
                   &           + zcoef1 * (  ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) )   &
-                  &                       - ( rhd(ji,jj,  jk)+rhd(ji,jj  ,jk-1) )  ) / e2v(ji,jj)
-!!gm                  &                       - ( rhd(ji,jj,  jk)+rhd(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj)
+                  &                       - ( rhd(ji,jj,  jk)+rhd(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj)
                ! add to the general momentum trend
                ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
@@ -305,8 +301,6 @@
             zcoef1 = zcoef0 * e3w_n(ji,jj,1)
             ! hydrostatic pressure gradient
-            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj  ,1) - rhd(ji,jj,1) ) / e1u(ji,jj)
-            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji  ,jj+1,1) - rhd(ji,jj,1) ) / e2v(ji,jj)
-!!gm            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj  ,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj)
-!!gm            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji  ,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj)
+            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj  ,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj)
+            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji  ,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj)
             ! add to the general momentum trend
             ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
@@ -324,11 +318,9 @@
                zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   &
                   &           + zcoef1 * (  ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) )   &
-                  &                       - ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) )  ) / e1u(ji,jj)
-!!gm                  &                       - ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) )  ) * r1_e1u(ji,jj)
+                  &                       - ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) )  ) * r1_e1u(ji,jj)
 
                zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)   &
                   &           + zcoef1 * (  ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) )   &
-                  &                       - ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) )  ) / e2v(ji,jj)
-!!gm                  &                       - ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj)
+                  &                       - ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj)
                ! add to the general momentum trend
                ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
@@ -349,6 +341,5 @@
                ua  (ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku)         ! subtract old value
                zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1)                   &   ! compute the new one
-                  &            + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) / e1u(ji,jj)
-!!gm                  &            + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) * r1_e1u(ji,jj)
+                  &            + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) * r1_e1u(ji,jj)
                ua  (ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku)         ! add the new one to the general momentum trend
             ENDIF
@@ -356,6 +347,5 @@
                va  (ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv)         ! subtract old value
                zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1)                   &   ! compute the new one
-                  &            + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) / e2v(ji,jj)
-!!gm                  &            + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) * r1_e2v(ji,jj)
+                  &            + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) * r1_e2v(ji,jj)
                va  (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv)         ! add the new one to the general momentum trend
             ENDIF
@@ -411,17 +401,13 @@
          DO ji = fs_2, fs_jpim1   ! vector opt.
             ! hydrostatic pressure gradient along s-surfaces
-!!gm            zhpi(ji,jj,1) = zcoef0 * r1_e1u(ji,jj) * ( e3w_n(ji+1,jj  ,1) * ( znad + rhd(ji+1,jj  ,1) )   &
-            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( e3w_n(ji+1,jj  ,1) * ( znad + rhd(ji+1,jj  ,1) )   &
-               &                                     - e3w_n(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) ) )
-!!gm            zhpj(ji,jj,1) = zcoef0 * r1_e2v(ji,jj) * ( e3w_n(ji  ,jj+1,1) * ( znad + rhd(ji  ,jj+1,1) )   &
-            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( e3w_n(ji  ,jj+1,1) * ( znad + rhd(ji  ,jj+1,1) )   &
-               &                                     - e3w_n(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) ) )
+            zhpi(ji,jj,1) = zcoef0 * r1_e1u(ji,jj) * (  e3w_n(ji+1,jj  ,1) * ( znad + rhd(ji+1,jj  ,1) )  &
+               &                                      - e3w_n(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) )  )
+            zhpj(ji,jj,1) = zcoef0 * r1_e2v(ji,jj) * (  e3w_n(ji  ,jj+1,1) * ( znad + rhd(ji  ,jj+1,1) )  &
+               &                                      - e3w_n(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) )  )
             ! s-coordinate pressure gradient correction
             zuap = -zcoef0 * ( rhd    (ji+1,jj,1) + rhd    (ji,jj,1) + 2._wp * znad )   &
-               &           * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) / e1u(ji,jj)
-!!gm               &           * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj)
+               &           * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj)
             zvap = -zcoef0 * ( rhd    (ji,jj+1,1) + rhd    (ji,jj,1) + 2._wp * znad )   &
-               &           * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) / e2v(ji,jj)
-!!gm               &           * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj)
+               &           * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj)
             ! add to the general momentum trend
             ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap
@@ -435,9 +421,7 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                ! hydrostatic pressure gradient along s-surfaces
-!!gm               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj)   &
-               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   &
+               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj)   &
                   &           * (  e3w_n(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )   &
                   &              - e3w_n(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad )  )
-!!gm               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   &
                zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj)   &
                   &           * (  e3w_n(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   &
@@ -445,9 +429,7 @@
                ! s-coordinate pressure gradient correction
                zuap = -zcoef0 * ( rhd    (ji+1,jj  ,jk) + rhd    (ji,jj,jk) + 2._wp * znad )   &
-                  &           * ( gde3w_n(ji+1,jj  ,jk) - gde3w_n(ji,jj,jk) ) / e1u(ji,jj)
-!!gm                  &           * ( gde3w_n(ji+1,jj  ,jk) - gde3w_n(ji,jj,jk) ) * r1_e1u(ji,jj)
+                  &           * ( gde3w_n(ji+1,jj  ,jk) - gde3w_n(ji,jj,jk) ) * r1_e1u(ji,jj)
                zvap = -zcoef0 * ( rhd    (ji  ,jj+1,jk) + rhd    (ji,jj,jk) + 2._wp * znad )   &
-                  &           * ( gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk) ) / e2v(ji,jj)
-!!gm                  &           * ( gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj)
+                  &           * ( gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj)
                ! add to the general momentum trend
                ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap
@@ -561,23 +543,19 @@
             ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure
             ! we assume ISF is in isostatic equilibrium
-!!gm            zhpi(ji,jj,1) = zcoef0 * r1_e1u(ji,jj) * ( 0.5_wp * e3w_n(ji+1,jj  ,iktp1i)                                    &
-            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * e3w_n(ji+1,jj  ,iktp1i)                                    &
-               &                                  * ( 2._wp * znad + rhd(ji+1,jj  ,iktp1i) + zrhdtop_oce(ji+1,jj  ) )   &
-               &                                  - 0.5_wp * e3w_n(ji  ,jj  ,ikt   )                                    &
-               &                                  * ( 2._wp * znad + rhd(ji  ,jj  ,ikt   ) + zrhdtop_oce(ji  ,jj  ) )   &
-               &                                  + ( ziceload(ji+1,jj) - ziceload(ji,jj))                              ) 
-!!gm            zhpj(ji,jj,1) = zcoef0 * r1_e2v(ji,jj) * ( 0.5_wp * e3w_n(ji  ,jj+1,iktp1j)                                    &
-            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w_n(ji  ,jj+1,iktp1j)                                    &
-               &                                  * ( 2._wp * znad + rhd(ji  ,jj+1,iktp1j) + zrhdtop_oce(ji  ,jj+1) )   &
-               &                                  - 0.5_wp * e3w_n(ji  ,jj  ,ikt   )                                    & 
-               &                                  * ( 2._wp * znad + rhd(ji  ,jj  ,ikt   ) + zrhdtop_oce(ji  ,jj  ) )   &
-               &                                  + ( ziceload(ji,jj+1) - ziceload(ji,jj) )                             ) 
+            zhpi(ji,jj,1) = zcoef0 * r1_e1u(ji,jj) * ( 0.5_wp * e3w_n(ji+1,jj  ,iktp1i)                                  &
+               &                                   * ( 2._wp * znad + rhd(ji+1,jj  ,iktp1i) + zrhdtop_oce(ji+1,jj  ) )   &
+               &                                   - 0.5_wp * e3w_n(ji  ,jj  ,ikt   )                                    &
+               &                                   * ( 2._wp * znad + rhd(ji  ,jj  ,ikt   ) + zrhdtop_oce(ji  ,jj  ) )   &
+               &                                   + ( ziceload(ji+1,jj) - ziceload(ji,jj) )                             ) 
+            zhpj(ji,jj,1) = zcoef0 * r1_e2v(ji,jj) * ( 0.5_wp * e3w_n(ji  ,jj+1,iktp1j)                                  &
+               &                                   * ( 2._wp * znad + rhd(ji  ,jj+1,iktp1j) + zrhdtop_oce(ji  ,jj+1) )   &
+               &                                   - 0.5_wp * e3w_n(ji  ,jj  ,ikt   )                                    & 
+               &                                   * ( 2._wp * znad + rhd(ji  ,jj  ,ikt   ) + zrhdtop_oce(ji  ,jj  ) )   &
+               &                                   + ( ziceload(ji,jj+1) - ziceload(ji,jj) )                             ) 
             ! s-coordinate pressure gradient correction (=0 if z coordinate)
             zuap = -zcoef0 * ( rhd    (ji+1,jj,1) + rhd    (ji,jj,1) + 2._wp * znad )   &
-               &           * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) / e1u(ji,jj)
-!!gm               &           * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj)
+               &           * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj)
             zvap = -zcoef0 * ( rhd    (ji,jj+1,1) + rhd    (ji,jj,1) + 2._wp * znad )   &
-               &           * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) / e2v(ji,jj)
-!!gm               &           * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj)
+               &           * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj)
             ! add to the general momentum trend
             ua(ji,jj,1) = ua(ji,jj,1) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1)
@@ -596,11 +574,9 @@
             IF ( iku .GT. 1 ) THEN
                ! case iku
-!!gm               zhpi(ji,jj,iku)   =  zcoef0 * r1_e1u(ji,jj) * ze3wu                                            &
-               zhpi(ji,jj,iku)   =  zcoef0 / e1u(ji,jj) * ze3wu                                            &
+               zhpi(ji,jj,iku)   =  zcoef0 * r1_e1u(ji,jj) * ze3wu                                         &
                       &                                 * ( rhd    (ji+1,jj,iku) + rhd   (ji,jj,iku)       &
                       &                                   + SIGN(1._wp,ze3wu) * grui(ji,jj) + 2._wp * znad )
                ! corrective term ( = 0 if z coordinate )
-!!gm               zuap              = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) * r1_e1u(ji,jj)
-               zuap              = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) / e1u(ji,jj)
+               zuap              = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) * r1_e1u(ji,jj)
                ! zhpi will be added in interior loop
                ua(ji,jj,iku)     = ua(ji,jj,iku) + zuap
@@ -609,12 +585,10 @@
 
                ! case iku + 1 (remove the zphi term added in the interior loop and compute the one corrected for zps)
-!!gm               zhpiint        =  zcoef0 * r1_e1u(ji,jj)                                                               &    
-               zhpiint        =  zcoef0 / e1u(ji,jj)                                                               &    
+               zhpiint        =  zcoef0 * r1_e1u(ji,jj)                                                               &    
                   &           * (  e3w_n(ji+1,jj  ,iku+1) * ( (rhd(ji+1,jj,iku+1) + znad)                          &
                   &                                         + (rhd(ji+1,jj,iku  ) + znad) ) * tmask(ji+1,jj,iku)   &
                   &              - e3w_n(ji  ,jj  ,iku+1) * ( (rhd(ji  ,jj,iku+1) + znad)                          &
                   &                                         + (rhd(ji  ,jj,iku  ) + znad) ) * tmask(ji  ,jj,iku)   )
-!!gm               zhpi(ji,jj,iku+1) =  zcoef0 * r1_e1u(ji,jj) * ge3rui(ji,jj) - zhpiint 
-               zhpi(ji,jj,iku+1) =  zcoef0 / e1u(ji,jj) * ge3rui(ji,jj) - zhpiint 
+               zhpi(ji,jj,iku+1) =  zcoef0 * r1_e1u(ji,jj) * ge3rui(ji,jj) - zhpiint 
             END IF
                
@@ -624,11 +598,9 @@
             IF ( ikv .GT. 1 ) THEN
                ! case ikv
-!!gm               zhpj(ji,jj,ikv)   =  zcoef0 * r1_e2v(ji,jj) * ze3wv                                            &
-               zhpj(ji,jj,ikv)   =  zcoef0 / e2v(ji,jj) * ze3wv                                            &
-                     &                                  * ( rhd(ji,jj+1,ikv) + rhd   (ji,jj,ikv)           &
-                     &                                    + SIGN(1._wp,ze3wv) * grvi(ji,jj) + 2._wp * znad )
+               zhpj(ji,jj,ikv)   =  zcoef0 * r1_e2v(ji,jj) * ze3wv                                            &
+                     &                                     * ( rhd(ji,jj+1,ikv) + rhd   (ji,jj,ikv)           &
+                     &                                       + SIGN(1._wp,ze3wv) * grvi(ji,jj) + 2._wp * znad )
                ! corrective term ( = 0 if z coordinate )
-!!gm               zvap              = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) * r1_e2v(ji,jj)
-               zvap              = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) / e2v(ji,jj)
+               zvap              = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) * r1_e2v(ji,jj)
                ! zhpi will be added in interior loop
                va(ji,jj,ikv)      = va(ji,jj,ikv) + zvap
@@ -637,12 +609,10 @@
                
                ! case ikv + 1 (remove the zphj term added in the interior loop and compute the one corrected for zps)
-!!gm               zhpjint        =  zcoef0 * r1_e2v(ji,jj)                                                              &
-               zhpjint        =  zcoef0 / e2v(ji,jj)                                                              &
+               zhpjint        =  zcoef0 * r1_e2v(ji,jj)                                                           &
                   &           * (  e3w_n(ji  ,jj+1,ikv+1) * ( (rhd(ji,jj+1,ikv+1) + znad)                         &
                   &                                       + (rhd(ji,jj+1,ikv  ) + znad) ) * tmask(ji,jj+1,ikv)    &
                   &              - e3w_n(ji  ,jj  ,ikv+1) * ( (rhd(ji,jj  ,ikv+1) + znad)                         &
                   &                                       + (rhd(ji,jj  ,ikv  ) + znad) ) * tmask(ji,jj  ,ikv)  )
-!!gm               zhpj(ji,jj,ikv+1) =  zcoef0 * r1_e2v(ji,jj) * ge3rvi(ji,jj) - zhpjint
-               zhpj(ji,jj,ikv+1) =  zcoef0 / e2v(ji,jj) * ge3rvi(ji,jj) - zhpjint
+               zhpj(ji,jj,ikv+1) =  zcoef0 * r1_e2v(ji,jj) * ge3rvi(ji,jj) - zhpjint
             END IF
          END DO
@@ -655,13 +625,9 @@
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
-!!gm useles !
-            iku=miku(ji,jj); ikv=mikv(ji,jj)
-!!gm
             DO jk = 2, jpkm1
                ! hydrostatic pressure gradient along s-surfaces
                ! zhpi is masked for the first wet cell  (contribution already done in the upper bloc)
                zhpi(ji,jj,jk) = zhpi(ji,jj,jk) + zhpi(ji,jj,jk-1)                                                     &
-!!gm                  &           + zcoef0 * r1_e1u(ji,jj)                                                                &
-                  &           + zcoef0 / e1u(ji,jj)                                                                &
+                  &           + zcoef0 * r1_e1u(ji,jj)                                                                &
                   &                    * ( e3w_n(ji+1,jj,jk) * ( (rhd(ji+1,jj,jk  ) + znad)                           &
                   &                                            + (rhd(ji+1,jj,jk-1) + znad) ) * tmask(ji+1,jj,jk-1)   &
@@ -671,6 +637,5 @@
                ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 
                zuap = - zcoef0 * ( rhd    (ji+1,jj  ,jk) + rhd    (ji,jj,jk) + 2._wp * znad )                    &
-                  &            * ( gde3w_n(ji+1,jj  ,jk) - gde3w_n(ji,jj,jk) ) / e1u(ji,jj) * umask(ji,jj,jk-1)
-!!gm                  &            * ( gde3w_n(ji+1,jj  ,jk) - gde3w_n(ji,jj,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk-1)
+                  &            * ( gde3w_n(ji+1,jj  ,jk) - gde3w_n(ji,jj,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk-1)
                ua(ji,jj,jk) = ua(ji,jj,jk) + ( zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk)
 
@@ -678,6 +643,5 @@
                ! zhpi is masked for the first wet cell  (contribution already done in the upper bloc)
                zhpj(ji,jj,jk) = zhpj(ji,jj,jk) + zhpj(ji,jj,jk-1)                                                     &
-!!gm                  &           + zcoef0 * r1_e2v(ji,jj)                                                                &
-                  &           + zcoef0 / e2v(ji,jj)                                                                &
+                  &           + zcoef0 * r1_e2v(ji,jj)                                                                &
                   &                    * ( e3w_n(ji  ,jj+1,jk) * ( (rhd(ji,jj+1,jk  ) + znad)                           &
                   &                                              + (rhd(ji,jj+1,jk-1) + znad) ) * tmask(ji,jj+1,jk-1)   &
@@ -687,6 +651,5 @@
                ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc)
                zvap = - zcoef0 * ( rhd    (ji  ,jj+1,jk) + rhd    (ji,jj,jk) + 2._wp * znad )                    &
-                  &            * ( gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk) ) / e2v(ji,jj) * vmask(ji,jj,jk-1)
-!!gm                  &            * ( gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk-1)
+                  &            * ( gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk-1)
                ! add to the general momentum trend
                va(ji,jj,jk) = va(ji,jj,jk) + ( zhpj(ji,jj,jk) + zvap ) * vmask(ji,jj,jk)
@@ -707,13 +670,10 @@
                ! remove old value (interior case)
                zuap            = -zcoef0 * ( rhd    (ji+1,jj  ,iku) + rhd    (ji,jj,iku) + 2._wp * znad )   &
-                     &                   * ( gde3w_n(ji+1,jj  ,iku) - gde3w_n(ji,jj,iku) ) / e1u(ji,jj)
-!!gm                     &                   * ( gde3w_n(ji+1,jj  ,iku) - gde3w_n(ji,jj,iku) ) * r1_e1u(ji,jj)
+                     &                   * ( gde3w_n(ji+1,jj  ,iku) - gde3w_n(ji,jj,iku) ) * r1_e1u(ji,jj)
                ua(ji,jj,iku)   = ua(ji,jj,iku) - zhpi(ji,jj,iku) - zuap
                ! put new value
                ! -zpshpi to avoid double contribution of the partial step in the top layer 
-!!gm               zuap            = -zcoef0 * ( aru(ji,jj) + 2._wp * znad ) * gzu(ji,jj)  * r1_e1u(ji,jj)
-               zuap            = -zcoef0 * ( aru(ji,jj) + 2._wp * znad ) * gzu(ji,jj)  / e1u(ji,jj)
-               zhpi(ji,jj,iku) =  zhpi(ji,jj,iku-1) + zcoef0 / e1u(ji,jj) * ge3ru(ji,jj) - zpshpi(ji,jj) 
-!!gm               zhpi(ji,jj,iku) =  zhpi(ji,jj,iku-1) + zcoef0 * r1_e1u(ji,jj) * ge3ru(ji,jj) - zpshpi(ji,jj) 
+               zuap            = -zcoef0 * ( aru(ji,jj) + 2._wp * znad ) * gzu(ji,jj)  * r1_e1u(ji,jj)
+               zhpi(ji,jj,iku) =  zhpi(ji,jj,iku-1) + zcoef0 * r1_e1u(ji,jj) * ge3ru(ji,jj) - zpshpi(ji,jj) 
                ua(ji,jj,iku)   =  ua(ji,jj,iku) + zhpi(ji,jj,iku) + zuap
             END IF
@@ -722,13 +682,10 @@
                ! remove old value (interior case)
                zvap            = -zcoef0 * ( rhd    (ji  ,jj+1,ikv) + rhd    (ji,jj,ikv) + 2._wp * znad )   &
-                     &                   * ( gde3w_n(ji  ,jj+1,ikv) - gde3w_n(ji,jj,ikv) )   / e2v(ji,jj)
-!!gm                     &                   * ( gde3w_n(ji  ,jj+1,ikv) - gde3w_n(ji,jj,ikv) )   * r1_e2v(ji,jj)
+                     &                   * ( gde3w_n(ji  ,jj+1,ikv) - gde3w_n(ji,jj,ikv) ) * r1_e2v(ji,jj)
                va(ji,jj,ikv)   = va(ji,jj,ikv) - zhpj(ji,jj,ikv) - zvap
                ! put new value
                ! -zpshpj to avoid double contribution of the partial step in the top layer 
-!!gm               zvap            = -zcoef0 * ( arv(ji,jj) + 2._wp * znad ) * gzv(ji,jj)     * r1_e2v(ji,jj)
-               zvap            = -zcoef0 * ( arv(ji,jj) + 2._wp * znad ) * gzv(ji,jj)     / e2v(ji,jj)
-               zhpj(ji,jj,ikv) =  zhpj(ji,jj,ikv-1) + zcoef0 / e2v(ji,jj) * ge3rv(ji,jj) - zpshpj(ji,jj)   
-!!gm               zhpj(ji,jj,ikv) =  zhpj(ji,jj,ikv-1) + zcoef0 * r1_e2v(ji,jj) * ge3rv(ji,jj) - zpshpj(ji,jj)   
+               zvap            = -zcoef0 * ( arv(ji,jj) + 2._wp * znad ) * gzv(ji,jj) * r1_e2v(ji,jj)
+               zhpj(ji,jj,ikv) =  zhpj(ji,jj,ikv-1) + zcoef0 * r1_e2v(ji,jj) * ge3rv(ji,jj) - zpshpj(ji,jj)   
                va(ji,jj,ikv)   =  va(ji,jj,ikv) + zhpj(ji,jj,ikv) + zvap
             END IF
@@ -938,8 +895,6 @@
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
-!!gm            zhpi(ji,jj,1) = ( rho_k(ji+1,jj  ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) * r1_e1u(ji,jj)
-            zhpi(ji,jj,1) = ( rho_k(ji+1,jj  ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) / e1u(ji,jj)
-            zhpj(ji,jj,1) = ( rho_k(ji  ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) / e2v(ji,jj)
-!!gm            zhpj(ji,jj,1) = ( rho_k(ji  ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) * r1_e2v(ji,jj)
+            zhpi(ji,jj,1) = ( rho_k(ji+1,jj  ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) * r1_e1u(ji,jj)
+            zhpj(ji,jj,1) = ( rho_k(ji  ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) * r1_e2v(ji,jj)
             ! add to the general momentum trend
             ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
@@ -957,10 +912,8 @@
                zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)                                &
                   &           + (  ( rho_k(ji+1,jj,jk) - rho_k(ji,jj,jk  ) )    &
-                  &              - ( rho_i(ji  ,jj,jk) - rho_i(ji,jj,jk-1) )  ) / e1u(ji,jj)
-!!gm                  &              - ( rho_i(ji  ,jj,jk) - rho_i(ji,jj,jk-1) )  ) * r1_e1u(ji,jj)
+                  &              - ( rho_i(ji  ,jj,jk) - rho_i(ji,jj,jk-1) )  ) * r1_e1u(ji,jj)
                zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)                                &
                   &           + (  ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk  ) )    &
-                  &               -( rho_j(ji,jj  ,jk) - rho_j(ji,jj,jk-1) )  ) / e2v(ji,jj)
-!!gm                  &               -( rho_j(ji,jj  ,jk) - rho_j(ji,jj,jk-1) )  ) * r1_e2v(ji,jj)
+                  &               -( rho_j(ji,jj  ,jk) - rho_j(ji,jj,jk-1) )  ) * r1_e2v(ji,jj)
                ! add to the general momentum trend
                ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
@@ -992,6 +945,6 @@
       !!
       INTEGER  ::   ji, jj, jk, jkk                 ! dummy loop indices
-      REAL(wp) ::   zcoef0, znad                    ! temporary scalars
-      !!
+      REAL(wp) ::   zcoef0, znad                    ! local scalars
+      !
       !! The local variables for the correction term
       INTEGER  :: jk1, jis, jid, jjs, jjd
@@ -1004,7 +957,7 @@
       !!----------------------------------------------------------------------
       !
-      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp )
-      CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh )
-      CALL wrk_alloc( jpi,jpj, zsshu_n, zsshv_n )
+      CALL wrk_alloc( jpi,jpj,jpk,   zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp )
+      CALL wrk_alloc( jpi,jpj,jpk,   zdept, zrhh )
+      CALL wrk_alloc( jpi,jpj,       zsshu_n, zsshv_n )
       !
       IF( kt == nit000 ) THEN
@@ -1014,5 +967,4 @@
       ENDIF
 
-      !!----------------------------------------------------------------------
       ! Local constant initialization
       zcoef0 = - grav
@@ -1060,12 +1012,11 @@
       ! 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)
+      CALL cspline( fsp, xsp, asp, bsp, csp, dsp, polynomial_type )
 
       ! Integrate the hydrostatic pressure "zhpi(:,:,:)" at "T(ji,jj,1)"
       DO jj = 2, jpj
         DO ji = 2, jpi
-          zrhdt1 = zrhh(ji,jj,1) - interp3(zdept(ji,jj,1),asp(ji,jj,1), &
-                                         bsp(ji,jj,1),   csp(ji,jj,1), &
-                                         dsp(ji,jj,1) ) * 0.25_wp * e3w_n(ji,jj,1)
+          zrhdt1 = zrhh(ji,jj,1) - interp3( zdept(ji,jj,1), asp(ji,jj,1), bsp(ji,jj,1),  &
+             &                                              csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w_n(ji,jj,1)
 
           ! assuming linear profile across the top half surface layer
@@ -1078,8 +1029,8 @@
         DO jj = 2, jpj
           DO ji = 2, jpi
-            zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) +                          &
-                             integ_spline(zdept(ji,jj,jk-1), zdept(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) +                                  &
+               &             integ_spline( zdept(ji,jj,jk-1), zdept(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
@@ -1091,4 +1042,10 @@
       DO jj = 2, jpjm1
         DO ji = 2, jpim1
+!!gm BUG ?    if it is ssh at u- & v-point then it should be:
+!          zsshu_n(ji,jj) = (e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji+1,jj) * sshn(ji+1,jj)) * &
+!                         & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 
+!          zsshv_n(ji,jj) = (e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji,jj+1) * sshn(ji,jj+1)) * &
+!                         & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 
+!!gm not this:
           zsshu_n(ji,jj) = (e1e2u(ji,jj) * sshn(ji,jj) + e1e2u(ji+1, jj) * sshn(ji+1,jj)) * &
                          & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 
@@ -1108,6 +1065,6 @@
         DO jj = 2, jpjm1
           DO ji = 2, jpim1
-            zu(ji,jj,jk) = zu(ji,jj,jk-1)- e3u_n(ji,jj,jk)
-            zv(ji,jj,jk) = zv(ji,jj,jk-1)- e3v_n(ji,jj,jk)
+            zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u_n(ji,jj,jk)
+            zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v_n(ji,jj,jk)
           END DO
         END DO
@@ -1126,8 +1083,8 @@
         DO jj = 2, jpjm1
           DO ji = 2, jpim1
-            zu(ji,jj,jk) = min(zu(ji,jj,jk), max(-zdept(ji,jj,jk), -zdept(ji+1,jj,jk)))
-            zu(ji,jj,jk) = max(zu(ji,jj,jk), min(-zdept(ji,jj,jk), -zdept(ji+1,jj,jk)))
-            zv(ji,jj,jk) = min(zv(ji,jj,jk), max(-zdept(ji,jj,jk), -zdept(ji,jj+1,jk)))
-            zv(ji,jj,jk) = max(zv(ji,jj,jk), min(-zdept(ji,jj,jk), -zdept(ji,jj+1,jk)))
+            zu(ji,jj,jk) = MIN(  zu(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) )  )
+            zu(ji,jj,jk) = MAX(  zu(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) )  )
+            zv(ji,jj,jk) = MIN(  zv(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) )  )
+            zv(ji,jj,jk) = MAX(  zv(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) )  )
           END DO
         END DO
@@ -1187,17 +1144,14 @@
                ! update the momentum trends in u direction
 
-!!gm               zdpdx1 = zcoef0 * r1_e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk))
-               zdpdx1 = zcoef0 / e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk))
+               zdpdx1 = zcoef0 * r1_e1u(ji,jj) * ( zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk) )
                IF( lk_vvl ) THEN
-!!gm                 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * &
-                 zdpdx2 = zcoef0 / e1u(ji,jj) * &
-                         ( REAL(jis-jid, wp) * (zpwes + zpwed) + (sshn(ji+1,jj)-sshn(ji,jj)) )
+                 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * &
+                    &    ( REAL(jis-jid, wp) * (zpwes + zpwed) + (sshn(ji+1,jj)-sshn(ji,jj)) )
                 ELSE
-!!gm                 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed)
-                 zdpdx2 = zcoef0 / e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed)
+                 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed)
                ENDIF
-
-               ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * &
-               &           umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk)
+!!gm  Since umask(ji,:,:) = tmask(ji,:,:) * tmask(ji+1,:,:)  by definition
+!!gm      in the line below only * umask(ji,jj,jk)  is needed !!
+               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,32 +1201,28 @@
                ! update the momentum trends in v direction
 
-!!gm               zdpdy1 = zcoef0 * r1_e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk))
-               zdpdy1 = zcoef0 / e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk))
+               zdpdy1 = zcoef0 * r1_e2v(ji,jj) * ( zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk) )
                IF( lk_vvl ) THEN
-!!gm                   zdpdy2 = zcoef0 * r1_e2v(ji,jj) * &
-                   zdpdy2 = zcoef0 / e2v(ji,jj) * &
+                   zdpdy2 = zcoef0 * r1_e2v(ji,jj) * &
                            ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (sshn(ji,jj+1)-sshn(ji,jj)) )
                ELSE
-!!gm                   zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd )
-                   zdpdy2 = zcoef0 / e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd )
+                   zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd )
                ENDIF
-
-               va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2)*&
-               &              vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk)
+!!gm  Since vmask(:,jj,:) = tmask(:,jj,:) * tmask(:,jj+1,:)  by definition
+!!gm      in the line below only * vmask(ji,jj,jk)  is needed !!
+               va(ji,jj,jk) = va(ji,jj,jk) + ( zdpdy1 + zdpdy2 ) * vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk)
             ENDIF
-
-
-           END DO
-        END DO
-      END DO
-      !
-      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp )
-      CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh )
-      CALL wrk_dealloc( jpi,jpj, zsshu_n, zsshv_n )
+               !
+            END DO
+         END DO
+      END DO
+      !
+      CALL wrk_dealloc( jpi,jpj,jpk,   zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp )
+      CALL wrk_dealloc( jpi,jpj,jpk,   zdept, zrhh )
+      CALL wrk_dealloc( jpi,jpj,       zsshu_n, zsshv_n )
       !
    END SUBROUTINE hpg_prj
 
 
-   SUBROUTINE cspline(fsp, xsp, asp, bsp, csp, dsp, polynomial_type)
+   SUBROUTINE cspline( fsp, xsp, asp, bsp, csp, dsp, polynomial_type )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE cspline  ***
@@ -1297,10 +1247,10 @@
       REAL(wp) ::   zdf(size(fsp,3))
       !!----------------------------------------------------------------------
-
+      !
+!!gm  WHAT !!!!!   THIS IS VERY DANGEROUS !!!!!  
       jpi   = size(fsp,1)
       jpj   = size(fsp,2)
       jpkm1 = size(fsp,3) - 1
-
-
+      !
       IF (polynomial_type == 1) THEN     ! Constrained Cubic Spline
          DO ji = 1, jpi
@@ -1335,8 +1285,7 @@
 
                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)
+                          &          ( 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)
+                          &          ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - 0.5_wp * zdf(jpkm1 - 1)
 
                DO jk = 1, jpkm1 - 1
@@ -1361,5 +1310,5 @@
          END DO
 
-      ELSE IF (polynomial_type == 2) THEN     ! Linear
+      ELSEIF ( polynomial_type == 2 ) THEN     ! Linear
          DO ji = 1, jpi
             DO jj = 1, jpj
@@ -1392,21 +1341,20 @@
       !!                extrapolation is also permitted (no value limit)
       !!----------------------------------------------------------------------
-      IMPLICIT NONE
       REAL(wp), INTENT(in) ::  x, xl, xr, fl, fr
       REAL(wp)             ::  f ! result of the interpolation (extrapolation)
       REAL(wp)             ::  zdeltx
       !!----------------------------------------------------------------------
-
+      !
       zdeltx = xr - xl
-      IF(abs(zdeltx) <= 10._wp * EPSILON(x)) THEN
-        f = 0.5_wp * (fl + fr)
+      IF( abs(zdeltx) <= 10._wp * EPSILON(x) ) THEN
+         f = 0.5_wp * (fl + fr)
       ELSE
-        f = ( (x - xl ) * fr - ( x - xr ) * fl ) / zdeltx
+         f = ( (x - xl ) * fr - ( x - xr ) * fl ) / zdeltx
       ENDIF
-
+      !
    END FUNCTION interp1
 
 
-   FUNCTION interp2(x, a, b, c, d)  RESULT(f)
+   FUNCTION interp2( x, a, b, c, d )  RESULT(f)
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE interp1  ***
@@ -1417,15 +1365,14 @@
       !!
       !!----------------------------------------------------------------------
-      IMPLICIT NONE
       REAL(wp), INTENT(in) ::  x, a, b, c, d
       REAL(wp)             ::  f ! value from the interpolation
       !!----------------------------------------------------------------------
-
+      !
       f = a + x* ( b + x * ( c + d * x ) )
-
+      !
    END FUNCTION interp2
 
 
-   FUNCTION interp3(x, a, b, c, d)  RESULT(f)
+   FUNCTION interp3( x, a, b, c, d )  RESULT(f)
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE interp1  ***
@@ -1437,15 +1384,14 @@
       !!
       !!----------------------------------------------------------------------
-      IMPLICIT NONE
       REAL(wp), INTENT(in) ::  x, a, b, c, d
       REAL(wp)             ::  f ! value from the interpolation
       !!----------------------------------------------------------------------
-
+      !
       f = b + x * ( 2._wp * c + 3._wp * d * x)
-
+      !
    END FUNCTION interp3
 
 
-   FUNCTION integ_spline(xl, xr, a, b, c, d)  RESULT(f)
+   FUNCTION integ_spline( xl, xr, a, b, c, d )  RESULT(f)
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE interp1  ***
@@ -1456,17 +1402,16 @@
       !!
       !!----------------------------------------------------------------------
-      IMPLICIT NONE
       REAL(wp), INTENT(in) ::  xl, xr, a, b, c, d
       REAL(wp)             ::  za1, za2, za3
       REAL(wp)             ::  f                   ! integration result
       !!----------------------------------------------------------------------
-
+      !
       za1 = 0.5_wp * b
       za2 = c / 3.0_wp
       za3 = 0.25_wp * d
-
+      !
       f  = xr * ( a + xr * ( za1 + xr * ( za2 + za3 * xr ) ) ) - &
          & xl * ( a + xl * ( za1 + xl * ( za2 + za3 * xl ) ) )
-
+      !
    END FUNCTION integ_spline
 
