Index: NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/tests/CANAL/MY_SRC/usrdef_istate.F90
===================================================================
--- NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/tests/CANAL/MY_SRC/usrdef_istate.F90	(revision 14037)
+++ NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/tests/CANAL/MY_SRC/usrdef_istate.F90	(revision 14062)
@@ -26,5 +26,6 @@
    PRIVATE
 
-   PUBLIC   usr_def_istate   ! called by istate.F90
+   PUBLIC   usr_def_istate       ! called by istate.F90
+   PUBLIC   usr_def_istate_ssh   ! called by sshwzv.F90
 
    !! * Substitutions
@@ -37,5 +38,5 @@
 CONTAINS
   
-   SUBROUTINE usr_def_istate( pdept, ptmask, pts, pu, pv, pssh )
+   SUBROUTINE usr_def_istate( pdept, ptmask, pts, pu, pv )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE usr_def_istate  ***
@@ -52,5 +53,4 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   pu      ! i-component of the velocity  [m/s] 
       REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   pv      ! j-component of the velocity  [m/s] 
-      REAL(wp), DIMENSION(jpi,jpj)         , INTENT(  out) ::   pssh    ! sea-surface height
       !
       INTEGER  :: ji, jj, jk, jl  ! dummy loop indices
@@ -87,7 +87,5 @@
 
       CASE(0)    ! rest
-         
-         ! sea level:
-         pssh(:,:) = 0.
+         !
          ! temperature:
          pts(:,:,:,jp_tem) = 10._wp
@@ -99,24 +97,5 @@
          
       CASE(1)    ! geostrophic zonal jet from -zjety to +zjety
-
-         ! sea level:
-         SELECT CASE( nn_fcase )
-         CASE(0)    ! f = f0
-            ! sea level: ssh = - fuy / g
-            WHERE( ABS(gphit) <= zjety )
-               pssh(:,:) = - ff_t(:,:) * rn_uzonal * gphit(:,:) * 1.e3 / grav
-            ELSEWHERE
-               pssh(:,:) = - ff_t(:,:) * rn_uzonal * SIGN(zjety, gphit(:,:)) * 1.e3 / grav
-            END WHERE
-         CASE(1)    ! f = f0 + beta*y
-            ! sea level: ssh = - u / g * ( fy + 0.5 * beta * y^2 )
-            zbeta = 2._wp * omega * COS( rad * rn_ppgphi0 ) / ra
-            WHERE( ABS(gphit) <= zjety )
-               pssh(:,:) = - rn_uzonal / grav * ( zf0 * gphit(:,:) * 1.e3 + 0.5 * zbeta * gphit(:,:) * gphit(:,:) * 1.e6 )
-            ELSEWHERE
-               pssh(:,:) = - rn_uzonal / grav * ( zf0 * SIGN(zjety, gphit(:,:)) * 1.e3   &
-                  &                             + 0.5 * zbeta * zjety * zjety * 1.e6 )
-            END WHERE
-         END SELECT
+         !
          ! temperature:
          pts(:,:,:,jp_tem) = 10._wp
@@ -139,25 +118,5 @@
          !                  
       CASE(2)    ! geostrophic zonal current shear
-      
-         ! sea level:
-         SELECT CASE( nn_fcase )
-         CASE(0)    ! f = f0
-            ! sea level: ssh = - fuy / g
-            WHERE( ABS(gphit) <= zjety )
-               pssh(:,:) = - ff_t(:,:) * rn_uzonal * ABS(gphit(:,:)) * 1.e3 / grav
-            ELSEWHERE
-               pssh(:,:) = - ff_t(:,:) * rn_uzonal * zjety * 1.e3 / grav
-            END WHERE
-         CASE(1)    ! f = f0 + beta*y
-            ! sea level: ssh = - u / g * ( fy + 0.5 * beta * y^2 )
-            zbeta = 2._wp * omega * COS( rad * rn_ppgphi0 ) / ra
-            WHERE( ABS(gphit) <= zjety )
-               pssh(:,:) = - SIGN(rn_uzonal, gphit(:,:)) / grav   &
-                  &        * ( zf0 * gphit(:,:) * 1.e3 + 0.5 * zbeta * gphit(:,:) * gphit(:,:) * 1.e6 )
-            ELSEWHERE
-               pssh(:,:) = - SIGN(rn_uzonal, gphit(:,:)) / grav   &
-                  &        * ( zf0 * SIGN(zjety, gphit(:,:)) * 1.e3 + 0.5 * zbeta * zjety * zjety * 1.e6 )
-            END WHERE
-         END SELECT
+         !
          ! temperature:
          pts(:,:,:,jp_tem) = 10._wp
@@ -176,5 +135,5 @@
          !                  
       CASE(3)    ! gaussian zonal currant
-
+         !
          ! zonal current
          DO jk=1, jpkm1
@@ -182,14 +141,4 @@
             pu(:,:,jk) = rn_uzonal * EXP( - 0.5 * gphit(:,:)**2 / rn_lambda**2 )
          END DO
-         
-         ! sea level:
-         pssh(:,1) = - ff_t(:,1) / grav * pu(:,1,1) * e2t(:,1)
-         DO jl=1, jpnj
-            DO_2D( 0, 0, 0, 0 )
-               pssh(ji,jj) = pssh(ji,jj-1) - ff_t(ji,jj) / grav * pu(ji,jj,1) * e2t(ji,jj)
-            END_2D
-            CALL lbc_lnk( 'usrdef_istate', pssh, 'T',  1. )
-         END DO
-         
          ! temperature:
          pts(:,:,:,jp_tem) = 10._wp
@@ -202,5 +151,5 @@
          !            
       CASE(4)    ! geostrophic zonal pulse
-   
+         !
          DO_2D( 1, 1, 1, 1 )
             IF ( ABS(glamt(ji,jj)) <= zjetx ) THEN
@@ -210,24 +159,22 @@
             ELSE
                zdu = 0.
-            END IF
+            ENDIF
             IF ( ABS(gphit(ji,jj)) <= zjety ) THEN
-               pssh(ji,jj) = - ff_t(ji,jj) * zdu * gphit(ji,jj) * 1.e3 / grav
                pu(ji,jj,:) = zdu
                pts(ji,jj,:,jp_sal) = zdu / rn_uzonal + 1.
             ELSE
-               pssh(ji,jj) = - ff_t(ji,jj) * zdu * SIGN(zjety,gphit(ji,jj)) * 1.e3 / grav 
                pu(ji,jj,:) = 0.
                pts(ji,jj,:,jp_sal) = 1.
-            END IF
-         END_2D
-         
+            ENDIF
+         END_2D
+         !
          ! temperature:
          pts(:,:,:,jp_tem) = 10._wp * ptmask(:,:,:)        
          pv(:,:,:) = 0.
-         
-       CASE(5)    ! vortex
-                  !
+         !
+      CASE(5)    ! vortex
+         !
          zf0   = 2._wp * omega * SIN( rad * rn_ppgphi0 )
-         zumax = rn_vtxmax * SIGN(1._wp, zf0) ! Here Anticyclonic: set zumax=-1 for cyclonic
+         zumax = rn_vtxmax * SIGN(1._wp, zf0)  ! Here Anticyclonic: set zumax=-1 for cyclonic
          zlambda = SQRT(2._wp)*rn_lambda*1.e3       ! Horizontal scale in meters 
          zn2 = 3.e-3**2
@@ -242,12 +189,4 @@
             ! Surface pressure: P(x,y,z) = F(z) * Psurf(x,y)
             zpsurf = zP0 * EXP(-(zx**2+zy**2)*zr_lambda2) - rho0 * ff_t(ji,jj) * rn_uzonal * zy
-            ! Sea level:
-            pssh(ji,jj) = 0.
-            DO jl=1,5
-               zdt = pssh(ji,jj)
-               zdzF = (1._wp - EXP(zdt-zH)) / (zH - 1._wp + EXP(-zH))   ! F'(z)
-               zrho1 = rho0 * (1._wp + zn2*zdt/grav) - zdzF * zpsurf / grav    ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y)
-               pssh(ji,jj) = zpsurf / (zrho1*grav) * ptmask(ji,jj,1)   ! ssh = Psurf / (Rho*g)
-            END DO
             ! temperature:
             DO jk=1,jpk
@@ -299,16 +238,146 @@
          !            
       END SELECT
-      
+      !
+      CALL lbc_lnk( 'usrdef_istate', pts , 'T',  1. )
+      CALL lbc_lnk_multi( 'usrdef_istate', pu, 'U', -1., pv, 'V', -1. )
+
+   END SUBROUTINE usr_def_istate
+
+  
+   SUBROUTINE usr_def_istate_ssh( ptmask, pssh )
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE usr_def_istate_ssh  ***
+      !! 
+      !! ** Purpose :   Initialization of the dynamics and tracers
+      !!                Here CANAL configuration 
+      !!
+      !! ** Method  :   Set ssh 
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   ptmask  ! t-point ocean mask             [m]
+      REAL(wp), DIMENSION(jpi,jpj)         , INTENT(  out) ::   pssh    ! sea-surface height
+      !
+      INTEGER  :: ji, jj, jk, jl  ! dummy loop indices
+      REAL(wp) :: zx, zy, zP0, zumax, zlambda, zr_lambda2, zn2, zf0, zH, zrho1, za, zf, zdzF
+      REAL(wp) :: zpsurf, zdyPs, zdxPs
+      REAL(wp) :: zdt, zdu, zdv
+      REAL(wp) :: zjetx, zjety, zbeta
+      REAL(wp), DIMENSION(jpi,jpj)  ::   zrandom
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'usr_def_istate_ssh : CANAL configuration, analytical definition of initial state'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   '
+      !
+      IF (ln_sshnoise) CALL RANDOM_NUMBER(zrandom)
+      zjetx = ABS(rn_ujetszx)/2.
+      zjety = ABS(rn_ujetszy)/2.
+      !
+      SELECT CASE(nn_initcase)
+      CASE(0)                      !==   rest  ==!
+         !
+         pssh(:,:) = 0.
+         !
+      CASE(1)                      !==  geostrophic zonal jet from -zjety to +zjety  ==!
+         !
+         SELECT CASE( nn_fcase )
+         CASE(0)                          !* f = f0 : ssh = - fuy / g
+            WHERE( ABS(gphit) <= zjety )
+               pssh(:,:) = - ff_t(:,:) * rn_uzonal * gphit(:,:) * 1.e3 / grav
+            ELSEWHERE
+               pssh(:,:) = - ff_t(:,:) * rn_uzonal * SIGN(zjety, gphit(:,:)) * 1.e3 / grav
+            END WHERE
+         CASE(1)                          !* f = f0 + beta*y : ssh = - u / g * ( fy + 0.5 * beta * y^2 )
+            zbeta = 2._wp * omega * COS( rad * rn_ppgphi0 ) / ra
+            WHERE( ABS(gphit) <= zjety )
+               pssh(:,:) = - rn_uzonal / grav * ( ff_t(:,:) * gphit(:,:) * 1.e3 + 0.5 * zbeta * gphit(:,:) * gphit(:,:) * 1.e6 )
+            ELSEWHERE
+               pssh(:,:) = - rn_uzonal / grav * ( ff_t(:,:) * SIGN(zjety, gphit(:,:)) * 1.e3   &
+                  &                             + 0.5 * zbeta * zjety * zjety * 1.e6 )
+            END WHERE
+         END SELECT
+         !                  
+      CASE(2)                      !==  geostrophic zonal current shear  ==!
+         !
+         SELECT CASE( nn_fcase )
+         CASE(0)                          !* f = f0 : ssh = - fuy / g
+            WHERE( ABS(gphit) <= zjety )
+               pssh(:,:) = - ff_t(:,:) * rn_uzonal * ABS(gphit(:,:)) * 1.e3 / grav
+            ELSEWHERE
+               pssh(:,:) = - ff_t(:,:) * rn_uzonal * zjety * 1.e3 / grav
+            END WHERE
+         CASE(1)                          !* f = f0 + beta*y : ssh = - u / g * ( fy + 0.5 * beta * y^2 )
+            zbeta = 2._wp * omega * COS( rad * rn_ppgphi0 ) / ra
+            WHERE( ABS(gphit) <= zjety )
+               pssh(:,:) = - SIGN(rn_uzonal, gphit(:,:)) / grav   &
+                  &        * ( ff_t(:,:) * gphit(:,:) * 1.e3 + 0.5 * zbeta * gphit(:,:) * gphit(:,:) * 1.e6 )
+            ELSEWHERE
+               pssh(:,:) = - SIGN(rn_uzonal, gphit(:,:)) / grav   &
+                  &        * ( ff_t(:,:) * SIGN(zjety, gphit(:,:)) * 1.e3 + 0.5 * zbeta * zjety * zjety * 1.e6 )
+            END WHERE
+         END SELECT
+         !                  
+      CASE(3)                      !==  gaussian zonal currant  ==!
+         !
+         pssh(:,1) = - ff_t(:,1) / grav * e2t(:,1) * rn_uzonal * EXP( - 0.5 * gphit(:,1)**2 / rn_lambda**2 )
+         DO jl=1, jpnj
+            DO_2D( 0, 0, 0, 0 )
+               pssh(ji,jj) = pssh(ji,jj-1) - ff_t(ji,jj) / grav * rn_uzonal * EXP( - 0.5 * gphit(ji,jj)**2 / rn_lambda**2 ) * e2t(ji,jj)
+            END_2D
+            CALL lbc_lnk( 'usrdef_istate_ssh', pssh, 'T',  1. )
+         END DO
+         !            
+      CASE(4)                      !==  geostrophic zonal pulse !!st need to implement a way to separate ssh properly  ==!
+         !
+         DO_2D( 1, 1, 1, 1 )
+            IF ( ABS(glamt(ji,jj)) <= zjetx ) THEN
+               zdu = rn_uzonal
+            ELSEIF ( ABS(glamt(ji,jj)) <= zjetx + 100. ) THEN
+               zdu = rn_uzonal * ( ( zjetx-ABS(glamt(ji,jj)) )/100. + 1. )
+            ELSE
+               zdu = 0.
+            ENDIF
+            IF ( ABS(gphit(ji,jj)) <= zjety ) THEN
+               pssh(ji,jj) = - ff_t(ji,jj) * zdu * gphit(ji,jj) * 1.e3 / grav
+            ELSE
+               pssh(ji,jj) = - ff_t(ji,jj) * zdu * SIGN(zjety,gphit(ji,jj)) * 1.e3 / grav 
+            ENDIF
+         END_2D
+         !
+      CASE(5)                    !==  vortex  ==!
+         !
+         zf0   = 2._wp * omega * SIN( rad * rn_ppgphi0 )
+         zumax = rn_vtxmax * SIGN(1._wp, zf0)   ! Here Anticyclonic: set zumax=-1 for cyclonic
+         zlambda = SQRT(2._wp)*rn_lambda        ! Horizontal scale in meters 
+         zn2 = 3.e-3**2
+         zH = 0.5_wp * 5000._wp
+         !
+         zr_lambda2 = 1._wp / zlambda**2
+         zP0 = rho0 * zf0 * zumax * zlambda * SQRT(EXP(1._wp)/2._wp)
+         !
+         DO_2D( 1, 1, 1, 1 )
+            zx = glamt(ji,jj) * 1.e3
+            zy = gphit(ji,jj) * 1.e3
+            !                                   ! Surface pressure: P(x,y,z) = F(z) * Psurf(x,y)
+            zpsurf = zP0 * EXP(-(zx**2+zy**2)*zr_lambda2) - rho0 * ff_t(ji,jj) * rn_uzonal * zy
+            pssh(ji,jj) = 0.
+            DO jl=1,5
+               zdt = pssh(ji,jj)
+               zdzF = (1._wp - EXP(zdt-zH)) / (zH - 1._wp + EXP(-zH))   ! F'(z)
+               zrho1 = rho0 * (1._wp + zn2*zdt/grav) - zdzF * zpsurf / grav    ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y)
+               pssh(ji,jj) = zpsurf / (zrho1*grav) * ptmask(ji,jj,1)   ! ssh = Psurf / (Rho*g)
+            END DO
+         END_2D
+         !            
+      END SELECT
+      !                          !==  add noise  ==!
       IF (ln_sshnoise) THEN
          CALL RANDOM_SEED()
          CALL RANDOM_NUMBER(zrandom)
          pssh(:,:) = pssh(:,:) + ( 0.1  * zrandom(:,:) - 0.05 )
-      END IF
-      CALL lbc_lnk( 'usrdef_istate', pssh, 'T',  1. )
-      CALL lbc_lnk( 'usrdef_istate', pts , 'T',  1. )
-      CALL lbc_lnk_multi( 'usrdef_istate', pu, 'U', -1., pv, 'V', -1. )
-
-   END SUBROUTINE usr_def_istate
-
+      ENDIF
+      CALL lbc_lnk( 'usrdef_istate_ssh', pssh, 'T',  1. )
+      !
+   END SUBROUTINE usr_def_istate_ssh
+   
    !!======================================================================
 END MODULE usrdef_istate
