Index: branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
===================================================================
--- branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 6484)
+++ branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 6579)
@@ -290,6 +290,5 @@
             !
             IF( ie1e2u_v == 0 ) THEN   ! e1e2u and e1e2v have not been read: compute them
-               !                          ! e2u and e1v does not include a
-               !                          reduction in some strait: apply reduction
+               !                       ! e2u and e1v does not include a reduction in some strait: apply reduction
                e1e2u (:,:) = e1u(:,:) * e2u(:,:)
                e1e2v (:,:) = e1v(:,:) * e2v(:,:)
@@ -301,5 +300,5 @@
             &                 gphit  , gphiu , gphiv , gphif ,    &
             &                 e1t    , e1u   , e1v   , e1f   ,    &
-            &                 e2t    , e2u   , e2v   , e2f   ,    & 
+            &                 e2t    , e2u   , e2v   , e2f   ,    &
             &                 ie1e2u_v  )
          !
@@ -331,24 +330,4 @@
       e2_e1u(:,:) = e2u(:,:) / e1u(:,:)
       e1_e2v(:,:) = e1v(:,:) / e2v(:,:)
-
-      IF( lwp .AND. nn_print >=1 .AND. .NOT.ln_rstart ) THEN      ! Control print : Grid informations (if not restart)
-         WRITE(numout,*)
-         WRITE(numout,*) '          longitude and e1 scale factors'
-         WRITE(numout,*) '          ------------------------------'
-         WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   &
-            glamv(ji,1), glamf(ji,1),   &
-            e1t(ji,1), e1u(ji,1),   &
-            e1v(ji,1), e1f(ji,1), ji = 1, jpi,10)
-9300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    &
-            f19.10, 1x, f19.10, 1x, f19.10, 1x, f19.10 )
-            !
-         WRITE(numout,*)
-         WRITE(numout,*) '          latitude and e2 scale factors'
-         WRITE(numout,*) '          -----------------------------'
-         WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   &
-            &                     gphiv(1,jj), gphif(1,jj),   &
-            &                     e2t  (1,jj), e2u  (1,jj),   &
-            &                     e2v  (1,jj), e2f  (1,jj), jj = 1, jpj, 10 )
-      ENDIF
 
 
Index: branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90
===================================================================
--- branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90	(revision 6484)
+++ branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90	(revision 6579)
@@ -71,13 +71,13 @@
       !
       !                                                         ! horizontal mesh (inum3)
-      CALL iom_rstput( 0, 0, inum0, 'glamt', glamt, ktype = jp_r4 )     !    ! latitude
-      CALL iom_rstput( 0, 0, inum0, 'glamu', glamu, ktype = jp_r4 )
-      CALL iom_rstput( 0, 0, inum0, 'glamv', glamv, ktype = jp_r4 )
-      CALL iom_rstput( 0, 0, inum0, 'glamf', glamf, ktype = jp_r4 )
-      
-      CALL iom_rstput( 0, 0, inum0, 'gphit', gphit, ktype = jp_r4 )     !    ! longitude
-      CALL iom_rstput( 0, 0, inum0, 'gphiu', gphiu, ktype = jp_r4 )
-      CALL iom_rstput( 0, 0, inum0, 'gphiv', gphiv, ktype = jp_r4 )
-      CALL iom_rstput( 0, 0, inum0, 'gphif', gphif, ktype = jp_r4 )
+      CALL iom_rstput( 0, 0, inum0, 'glamt', glamt, ktype = jp_r8 )     !    ! latitude
+      CALL iom_rstput( 0, 0, inum0, 'glamu', glamu, ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum0, 'glamv', glamv, ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum0, 'glamf', glamf, ktype = jp_r8 )
+      
+      CALL iom_rstput( 0, 0, inum0, 'gphit', gphit, ktype = jp_r8 )     !    ! longitude
+      CALL iom_rstput( 0, 0, inum0, 'gphiu', gphiu, ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum0, 'gphiv', gphiv, ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum0, 'gphif', gphif, ktype = jp_r8 )
       
       CALL iom_rstput( 0, 0, inum0, 'e1t', e1t, ktype = jp_r8 )         !    ! e1 scale factors
@@ -229,13 +229,13 @@
 
       !                                                         ! horizontal mesh (inum3)
-      CALL iom_rstput( 0, 0, inum3, 'glamt', glamt, ktype = jp_r4 )     !    ! latitude
-      CALL iom_rstput( 0, 0, inum3, 'glamu', glamu, ktype = jp_r4 )
-      CALL iom_rstput( 0, 0, inum3, 'glamv', glamv, ktype = jp_r4 )
-      CALL iom_rstput( 0, 0, inum3, 'glamf', glamf, ktype = jp_r4 )
-      
-      CALL iom_rstput( 0, 0, inum3, 'gphit', gphit, ktype = jp_r4 )     !    ! longitude
-      CALL iom_rstput( 0, 0, inum3, 'gphiu', gphiu, ktype = jp_r4 )
-      CALL iom_rstput( 0, 0, inum3, 'gphiv', gphiv, ktype = jp_r4 )
-      CALL iom_rstput( 0, 0, inum3, 'gphif', gphif, ktype = jp_r4 )
+      CALL iom_rstput( 0, 0, inum3, 'glamt', glamt, ktype = jp_r8 )     !    ! latitude
+      CALL iom_rstput( 0, 0, inum3, 'glamu', glamu, ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum3, 'glamv', glamv, ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum3, 'glamf', glamf, ktype = jp_r8 )
+      
+      CALL iom_rstput( 0, 0, inum3, 'gphit', gphit, ktype = jp_r8 )     !    ! longitude
+      CALL iom_rstput( 0, 0, inum3, 'gphiu', gphiu, ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum3, 'gphiv', gphiv, ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum3, 'gphif', gphif, ktype = jp_r8 )
       
       CALL iom_rstput( 0, 0, inum3, 'e1t', e1t, ktype = jp_r8 )         !    ! e1 scale factors
@@ -253,9 +253,9 @@
       ! note that mbkt is set to 1 over land ==> use surface tmask
       zprt(:,:) = ssmask(:,:) * REAL( mbkt(:,:) , wp )
-      CALL iom_rstput( 0, 0, inum4, 'mbathy', zprt, ktype = jp_i2 )     !    ! nb of ocean T-points
+      CALL iom_rstput( 0, 0, inum4, 'mbathy', zprt, ktype = jp_i4 )     !    ! nb of ocean T-points
       zprt(:,:) = ssmask(:,:) * REAL( mikt(:,:) , wp )
-      CALL iom_rstput( 0, 0, inum4, 'misf', zprt, ktype = jp_i2 )       !    ! nb of ocean T-points
+      CALL iom_rstput( 0, 0, inum4, 'misf', zprt, ktype = jp_i4 )       !    ! nb of ocean T-points
       zprt(:,:) = ssmask(:,:) * REAL( risfdep(:,:) , wp )
-      CALL iom_rstput( 0, 0, inum4, 'isfdraft', zprt, ktype = jp_r4 )       !    ! nb of ocean T-points
+      CALL iom_rstput( 0, 0, inum4, 'isfdraft', zprt, ktype = jp_r8 )       !    ! nb of ocean T-points
             
       IF( ln_sco ) THEN                                         ! s-coordinate
@@ -279,6 +279,6 @@
          CALL iom_rstput( 0, 0, inum4, 'gdept_1d' , gdept_1d )  !    ! stretched system
          CALL iom_rstput( 0, 0, inum4, 'gdepw_1d' , gdepw_1d )
-         CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0, ktype = jp_r4 )     
-         CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0, ktype = jp_r4 )
+         CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0, ktype = jp_r8 )
+         CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0, ktype = jp_r8 )
       ENDIF
       
@@ -302,5 +302,5 @@
          !
          IF( nmsh <= 3 ) THEN                                   !    ! 3D depth
-            CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0, ktype = jp_r4 )     
+            CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0, ktype = jp_r8 )
             DO jk = 1,jpk   
                DO jj = 1, jpjm1   
@@ -308,11 +308,11 @@
                      zdepu(ji,jj,jk) = MIN( gdept_0(ji,jj,jk) , gdept_0(ji+1,jj  ,jk) )
                      zdepv(ji,jj,jk) = MIN( gdept_0(ji,jj,jk) , gdept_0(ji  ,jj+1,jk) )
-                  END DO   
-               END DO   
+                  END DO
+               END DO
             END DO
-            CALL lbc_lnk( zdepu, 'U', 1. )   ;   CALL lbc_lnk( zdepv, 'V', 1. ) 
-            CALL iom_rstput( 0, 0, inum4, 'gdepu', zdepu, ktype = jp_r4 )
-            CALL iom_rstput( 0, 0, inum4, 'gdepv', zdepv, ktype = jp_r4 )
-            CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0, ktype = jp_r4 )
+            CALL lbc_lnk( zdepu, 'U', 1. )   ;   CALL lbc_lnk( zdepv, 'V', 1. )
+            CALL iom_rstput( 0, 0, inum4, 'gdepu', zdepu, ktype = jp_r8 )
+            CALL iom_rstput( 0, 0, inum4, 'gdepv', zdepv, ktype = jp_r8 )
+            CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0, ktype = jp_r8 )
          ELSE                                                   !    ! 2D bottom depth
             DO jj = 1,jpj   
@@ -322,6 +322,6 @@
                END DO
             END DO
-            CALL iom_rstput( 0, 0, inum4, 'hdept', zprt, ktype = jp_r4 )     
-            CALL iom_rstput( 0, 0, inum4, 'hdepw', zprw, ktype = jp_r4 ) 
+            CALL iom_rstput( 0, 0, inum4, 'hdept', zprt, ktype = jp_r8 )
+            CALL iom_rstput( 0, 0, inum4, 'hdepw', zprw, ktype = jp_r8 )
          ENDIF
          !
Index: branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90
===================================================================
--- branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90	(revision 6484)
+++ branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90	(revision 6579)
@@ -102,4 +102,5 @@
       IF(  ln_tsd_init .OR. ln_tsd_tradmp  ) THEN
          !
+IF(lwp) WRITE(numout,*) ' sto dentro if ln_tsd_init prima di allocate if necessary riga 104 '
          ALLOCATE( sf_tsd(jpts), STAT=ierr0 )
          IF( ierr0 > 0 ) THEN
@@ -107,6 +108,8 @@
          ENDIF
          !
+IF(lwp) WRITE(numout,*) ' riga 110 : allocate jp_temp jpk '
                                 ALLOCATE( sf_tsd(jp_tem)%fnow(jpi,jpj,jpk)   , STAT=ierr0 )
          IF( sn_tem%ln_tint )   ALLOCATE( sf_tsd(jp_tem)%fdta(jpi,jpj,jpk,2) , STAT=ierr1 )
+IF(lwp) WRITE(numout,*) ' riga 113 : dopo allocate jp_temp jpk, 2 '
                                 ALLOCATE( sf_tsd(jp_sal)%fnow(jpi,jpj,jpk)   , STAT=ierr2 )
          IF( sn_sal%ln_tint )   ALLOCATE( sf_tsd(jp_sal)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 )
@@ -117,5 +120,7 @@
          !                         ! fill sf_tsd with sn_tem & sn_sal and control print
          slf_i(jp_tem) = sn_tem   ;   slf_i(jp_sal) = sn_sal
+IF(lwp) WRITE(numout,*) ' riga 122 : prima di fld_fill'
          CALL fld_fill( sf_tsd, slf_i, cn_dir, 'dta_tsd', 'Temperature & Salinity data', 'namtsd' )
+IF(lwp) WRITE(numout,*) ' riga 124 : dopo di fld_fill'
          !
       ENDIF
Index: branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
===================================================================
--- branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 6484)
+++ branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 6579)
@@ -14,12 +14,10 @@
    !!            3.3  !  2010-10  (C. Ethe) merge TRC-TRA
    !!            3.4  !  2011-04  (G. Madec) Merge of dtatem and dtasal & suppression of tb,tn/sb,sn 
+   !!            3.7  !  2016-04  (S. Flavoni) change configuration's interface:
+   !!                              read file or CALL usr_def module to compute initial state (example given for GYRE)
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
    !!   istate_init   : initial state setting
-   !!   istate_tem    : analytical profile for initial Temperature
-   !!   istate_sal    : analytical profile for initial Salinity
-   !!   istate_eel    : initial state setting of EEL R5 configuration
-   !!   istate_gyre   : initial state setting of GYRE configuration
    !!   istate_uvg    : initial velocity in geostropic balance
    !!----------------------------------------------------------------------
@@ -43,4 +41,7 @@
    USE wrk_nemo        ! Memory allocation
    USE timing          ! Timing
+   USE usrdef          ! User defined routine
+
+
 
    IMPLICIT NONE
@@ -48,4 +49,5 @@
 
    PUBLIC   istate_init   ! routine called by step.F90
+   !SF PUBLIC   ini_read      ! subroutine ini_read 
 
    !! * Substitutions
@@ -75,4 +77,5 @@
       IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
 
+     !SF initialisation of T & S with data file
                      CALL dta_tsd_init        ! Initialisation of T & S input data
       IF( lk_c1d )   CALL dta_uvd_init        ! Initialization of U & V input data
@@ -91,8 +94,8 @@
          !                                    ! Start from rest
          !                                    ! ---------------
-         numror = 0                              ! define numror = 0 -> no restart file to read
-         neuler = 0                              ! Set time-step indicator at nit000 (euler forward)
-         CALL day_init                           ! model calendar (using both namelist and restart infos)
-         !                                       ! Initialization of ocean to zero
+         numror = 0                           ! define numror = 0 -> no restart file to read
+         neuler = 0                           ! Set time-step indicator at nit000 (euler forward)
+         CALL day_init                        ! model calendar (using both namelist and restart infos)
+         !                                    ! Initialization of ocean to zero
          !   before fields      !       now fields     
          sshb (:,:)   = 0._wp   ;   sshn (:,:)   = 0._wp
@@ -101,23 +104,17 @@
                                     hdivn(:,:,:) = 0._wp
          !
-         IF( cp_cfg == 'eel' ) THEN
-            CALL istate_eel                      ! EEL   configuration : start from pre-defined U,V T-S fields
-         ELSEIF( cp_cfg == 'gyre' ) THEN         
-            CALL istate_gyre                     ! GYRE  configuration : start from pre-defined T-S fields
-         ELSE                                    ! Initial T-S, U-V fields read in files
-            IF ( ln_tsd_init ) THEN              ! read 3D T and S data at nit000
-               CALL dta_tsd( nit000, tsb )  
-               tsn(:,:,:,:) = tsb(:,:,:,:)
-               !
-            ELSE                                 ! Initial T-S fields defined analytically
-               CALL istate_t_s
-            ENDIF
-            IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000
-               CALL wrk_alloc( jpi,jpj,jpk,2,   zuvd )
-               CALL dta_uvd( nit000, zuvd )
-               ub(:,:,:) = zuvd(:,:,:,1) ;  un(:,:,:) = ub(:,:,:)
-               vb(:,:,:) = zuvd(:,:,:,2) ;  vn(:,:,:) = vb(:,:,:)
-               CALL wrk_dealloc( jpi,jpj,jpk,2,   zuvd )
-            ENDIF
+         IF( ln_tsd_init ) THEN               ! read 3D T and S data at nit000
+            CALL dta_tsd( nit000, tsb )  
+         ELSE                                 ! user defined initial T and S
+            CALL usr_def_ini( tsb )                
+         ENDIF
+         tsn(:,:,:,:) = tsb(:,:,:,:)          ! set now to before values
+         !
+         IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000
+            CALL wrk_alloc( jpi,jpj,jpk,2,   zuvd )
+            CALL dta_uvd( nit000, zuvd )
+            ub(:,:,:) = zuvd(:,:,:,1) ;  un(:,:,:) = ub(:,:,:)
+            vb(:,:,:) = zuvd(:,:,:,2) ;  vn(:,:,:) = vb(:,:,:)
+            CALL wrk_dealloc( jpi,jpj,jpk,2,   zuvd )
          ENDIF
          !
@@ -131,5 +128,5 @@
 !!gm 
          ! 
-      ENDIF
+      ENDIF 
       ! 
       ! Initialize "now" and "before" barotropic velocities:
@@ -162,345 +159,3 @@
    END SUBROUTINE istate_init
 
-
-   SUBROUTINE istate_t_s
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE istate_t_s  ***
-      !!   
-      !! ** Purpose :   Intialization of the temperature field with an 
-      !!      analytical profile or a file (i.e. in EEL configuration)
-      !!
-      !! ** Method  : - temperature: use Philander analytic profile
-      !!              - salinity   : use to a constant value 35.5
-      !!
-      !! References :  Philander ???
-      !!----------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   zsal = 35.50_wp
-      !!----------------------------------------------------------------------
-      !
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'istate_t_s : Philander s initial temperature profile'
-      IF(lwp) WRITE(numout,*) '~~~~~~~~~~   and constant salinity (',zsal,' psu)'
-      !
-      DO jk = 1, jpk
-         tsn(:,:,jk,jp_tem) = (  ( ( 7.5 - 0. * ABS( gphit(:,:) )/30. ) * ( 1.-TANH((gdept_n(:,:,jk)-80.)/30.) )   &
-            &                + 10. * ( 5000. - gdept_n(:,:,jk) ) /5000.)  ) * tmask(:,:,jk)
-         tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem)
-      END DO
-      tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:)
-      tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
-      !
-   END SUBROUTINE istate_t_s
-
-
-   SUBROUTINE istate_eel
-      !!----------------------------------------------------------------------
-      !!                   ***  ROUTINE istate_eel  ***
-      !! 
-      !! ** Purpose :   Initialization of the dynamics and tracers for EEL R5
-      !!      configuration (channel with or without a topographic bump)
-      !!
-      !! ** Method  : - set temprature field
-      !!              - set salinity field
-      !!              - set velocity field including horizontal divergence
-      !!                and relative vorticity fields
-      !!----------------------------------------------------------------------
-      USE divhor     ! hor. divergence      (div_hor routine)
-      USE iom
-      !
-      INTEGER  ::   inum              ! temporary logical unit
-      INTEGER  ::   ji, jj, jk        ! dummy loop indices
-      INTEGER  ::   ijloc
-      REAL(wp) ::   zh1, zh2, zslope, zcst, zfcor   ! temporary scalars
-      REAL(wp) ::   zt1  = 15._wp                   ! surface temperature value (EEL R5)
-      REAL(wp) ::   zt2  =  5._wp                   ! bottom  temperature value (EEL R5)
-      REAL(wp) ::   zsal = 35.0_wp                  ! constant salinity (EEL R2, R5 and R6)
-      REAL(wp) ::   zueel = 0.1_wp                  ! constant uniform zonal velocity (EEL R5)
-      REAL(wp), DIMENSION(jpiglo,jpjglo) ::   zssh  ! initial ssh over the global domain
-      !!----------------------------------------------------------------------
-      !
-      SELECT CASE ( jp_cfg ) 
-         !                                              ! ====================
-         CASE ( 5 )                                     ! EEL R5 configuration
-            !                                           ! ====================
-            !
-            ! set temperature field with a linear profile
-            ! -------------------------------------------
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: linear temperature profile'
-            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
-            !
-            zh1 = gdept_1d(  1  )
-            zh2 = gdept_1d(jpkm1)
-            !
-            zslope = ( zt1 - zt2 ) / ( zh1 - zh2 )
-            zcst   = ( zt1 * ( zh1 - zh2) - ( zt1 - zt2 ) * zh1 ) / ( zh1 - zh2 )
-            !
-            DO jk = 1, jpk
-               tsn(:,:,jk,jp_tem) = ( zt2 + zt1 * exp( - gdept_n(:,:,jk) / 1000 ) ) * tmask(:,:,jk)
-               tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem)
-            END DO
-            !
-            ! set salinity field to a constant value
-            ! --------------------------------------
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: constant salinity field, S = ', zsal
-            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
-            !
-            tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:)
-            tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
-            !
-            ! set the dynamics: U,V, hdiv (and ssh if necessary)
-            ! ----------------
-            ! Start EEL5 configuration with barotropic geostrophic velocities 
-            ! according the sshb and sshn SSH imposed.
-            ! we assume a uniform grid (hence the use of e1t(1,1) for delta_y)
-            ! we use the Coriolis frequency at mid-channel.   
-            ub(:,:,:) = zueel * umask(:,:,:)
-            un(:,:,:) = ub(:,:,:)
-            ijloc = mj0(INT(jpjglo-1)/2)
-            zfcor = ff(1,ijloc)
-            !
-            DO jj = 1, jpjglo
-               zssh(:,jj) = - (FLOAT(jj)- FLOAT(jpjglo-1)/2.)*zueel*e1t(1,1)*zfcor/grav 
-            END DO
-            !
-            IF(lwp) THEN
-               WRITE(numout,*) ' Uniform zonal velocity for EEL R5:',zueel
-               WRITE(numout,*) ' Geostrophic SSH profile as a function of y:'
-               WRITE(numout,'(12(1x,f6.2))') zssh(1,:)
-            ENDIF
-            !
-            DO jj = 1, nlcj
-               DO ji = 1, nlci
-                  sshb(ji,jj) = zssh( mig(ji) , mjg(jj) ) * tmask(ji,jj,1)
-               END DO
-            END DO
-            sshb(nlci+1:jpi,      :   ) = 0.e0      ! set to zero extra mpp columns
-            sshb(      :   ,nlcj+1:jpj) = 0.e0      ! set to zero extra mpp rows
-            !
-            sshn(:,:) = sshb(:,:)                   ! set now ssh to the before value
-            !
-            IF( nn_rstssh /= 0 ) THEN  
-               nn_rstssh = 0                        ! hand-made initilization of ssh 
-               CALL ctl_warn( 'istate_eel: force nn_rstssh = 0' )
-            ENDIF
-            !
-!!gm  Check  here call to div_hor should not be necessary
-!!gm         div_hor call runoffs  not sure they are defined at that level
-            CALL div_hor( nit000 )                  ! horizontal divergence and relative vorticity (curl)
-            ! N.B. the vertical velocity will be computed from the horizontal divergence field
-            ! in istate by a call to wzv routine
-
-
-            !                                     ! ==========================
-         CASE ( 2 , 6 )                           ! EEL R2 or R6 configuration
-            !                                     ! ==========================
-            !
-            ! set temperature field with a NetCDF file
-            ! ----------------------------------------
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) 'istate_eel : EEL R2 or R6: read initial temperature in a NetCDF file'
-            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
-            !
-            CALL iom_open ( 'eel.initemp', inum )
-            CALL iom_get ( inum, jpdom_data, 'initemp', tsb(:,:,:,jp_tem) ) ! read before temprature (tb)
-            CALL iom_close( inum )
-            !
-            tsn(:,:,:,jp_tem) = tsb(:,:,:,jp_tem)                            ! set nox temperature to tb
-            !
-            ! set salinity field to a constant value
-            ! --------------------------------------
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: constant salinity field, S = ', zsal
-            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
-            !
-            tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:)
-            tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
-            !
-            !                                    ! ===========================
-         CASE DEFAULT                            ! NONE existing configuration
-            !                                    ! ===========================
-            WRITE(ctmp1,*) 'EEL with a ', jp_cfg,' km resolution is not coded'
-            CALL ctl_stop( ctmp1 )
-            !
-      END SELECT
-      !
-   END SUBROUTINE istate_eel
-
-
-   SUBROUTINE istate_gyre
-      !!----------------------------------------------------------------------
-      !!                   ***  ROUTINE istate_gyre  ***
-      !! 
-      !! ** Purpose :   Initialization of the dynamics and tracers for GYRE
-      !!      configuration (double gyre with rotated domain)
-      !!
-      !! ** Method  : - set temprature field
-      !!              - set salinity   field
-      !!----------------------------------------------------------------------
-      INTEGER :: ji, jj, jk  ! dummy loop indices
-      INTEGER            ::   inum          ! temporary logical unit
-      INTEGER, PARAMETER ::   ntsinit = 0   ! (0/1) (analytical/input data files) T&S initialization
-      !!----------------------------------------------------------------------
-      !
-      SELECT CASE ( ntsinit)
-      !
-      CASE ( 0 )                  ! analytical T/S profil deduced from LEVITUS
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'istate_gyre : initial analytical T and S profil deduced from LEVITUS '
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
-         !
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  tsn(ji,jj,jk,jp_tem) = (  16. - 12. * TANH( (gdept_n(ji,jj,jk) - 400) / 700 )         )   &
-                       &           * (-TANH( (500-gdept_n(ji,jj,jk)) / 150 ) + 1) / 2               &
-                       &       + (      15. * ( 1. - TANH( (gdept_n(ji,jj,jk)-50.) / 1500.) )       &
-                       &                - 1.4 * TANH((gdept_n(ji,jj,jk)-100.) / 100.)               &    
-                       &                + 7.  * (1500. - gdept_n(ji,jj,jk)) / 1500.             )   & 
-                       &           * (-TANH( (gdept_n(ji,jj,jk) - 500) / 150) + 1) / 2
-                  tsn(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
-                  tsb(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem)
-
-                  tsn(ji,jj,jk,jp_sal) =  (  36.25 - 1.13 * TANH( (gdept_n(ji,jj,jk) - 305) / 460 )  )  &
-                     &              * (-TANH((500 - gdept_n(ji,jj,jk)) / 150) + 1) / 2          &
-                     &          + (  35.55 + 1.25 * (5000. - gdept_n(ji,jj,jk)) / 5000.         &
-                     &                - 1.62 * TANH( (gdept_n(ji,jj,jk) - 60.  ) / 650. )       &
-                     &                + 0.2  * TANH( (gdept_n(ji,jj,jk) - 35.  ) / 100. )       &
-                     &                + 0.2  * TANH( (gdept_n(ji,jj,jk) - 1000.) / 5000.)    )  &
-                     &              * (-TANH((gdept_n(ji,jj,jk) - 500) / 150) + 1) / 2 
-                  tsn(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
-                  tsb(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal)
-               END DO
-            END DO
-         END DO
-         !
-      CASE ( 1 )                  ! T/S data fields read in dta_tem.nc/data_sal.nc files
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'istate_gyre : initial T and S read from dta_tem.nc/data_sal.nc files'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
-         IF(lwp) WRITE(numout,*) '              NetCDF FORMAT'
-
-         ! Read temperature field
-         ! ----------------------
-         CALL iom_open ( 'data_tem', inum )
-         CALL iom_get ( inum, jpdom_data, 'votemper', tsn(:,:,:,jp_tem) ) 
-         CALL iom_close( inum )
-
-         tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:) 
-         tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)
-
-         ! Read salinity field
-         ! -------------------
-         CALL iom_open ( 'data_sal', inum )
-         CALL iom_get ( inum, jpdom_data, 'vosaline', tsn(:,:,:,jp_sal) ) 
-         CALL iom_close( inum )
-
-         tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:) 
-         tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
-         !
-      END SELECT
-      !
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) '              Initial temperature and salinity profiles:'
-         WRITE(numout, "(9x,' level   gdept_1d   temperature   salinity   ')" )
-         WRITE(numout, "(10x, i4, 3f10.2)" ) ( jk, gdept_1d(jk), tsn(2,2,jk,jp_tem), tsn(2,2,jk,jp_sal), jk = 1, jpk )
-      ENDIF
-      !
-   END SUBROUTINE istate_gyre
-
-
-   SUBROUTINE istate_uvg
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE istate_uvg  ***
-      !!
-      !! ** Purpose :   Compute the geostrophic velocities from (tn,sn) fields
-      !!
-      !! ** Method  :   Using the hydrostatic hypothesis the now hydrostatic 
-      !!      pressure is computed by integrating the in-situ density from the
-      !!      surface to the bottom.
-      !!                 p=integral [ rau*g dz ]
-      !!----------------------------------------------------------------------
-      USE divhor          ! hor. divergence                       (div_hor routine)
-      USE lbclnk          ! ocean lateral boundary condition (or mpp link)
-      !
-      INTEGER ::   ji, jj, jk        ! dummy loop indices
-      REAL(wp) ::   zmsv, zphv, zmsu, zphu, zalfg     ! temporary scalars
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zprn
-      !!----------------------------------------------------------------------
-      !
-      CALL wrk_alloc( jpi,jpj,jpk,   zprn)
-      !
-      IF(lwp) WRITE(numout,*) 
-      IF(lwp) WRITE(numout,*) 'istate_uvg : Start from Geostrophy'
-      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
-
-      ! Compute the now hydrostatic pressure
-      ! ------------------------------------
-
-      zalfg = 0.5 * grav * rau0
-      
-      zprn(:,:,1) = zalfg * e3w_n(:,:,1) * ( 1 + rhd(:,:,1) )       ! Surface value
-
-      DO jk = 2, jpkm1                                              ! Vertical integration from the surface
-         zprn(:,:,jk) = zprn(:,:,jk-1)   &
-            &         + zalfg * e3w_n(:,:,jk) * ( 2. + rhd(:,:,jk) + rhd(:,:,jk-1) )
-      END DO  
-
-      ! Compute geostrophic balance
-      ! ---------------------------
-      DO jk = 1, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vertor opt.
-               zmsv = 1. / MAX(  umask(ji-1,jj+1,jk) + umask(ji  ,jj+1,jk)   &
-                               + umask(ji-1,jj  ,jk) + umask(ji  ,jj  ,jk) , 1.  )
-               zphv = ( zprn(ji  ,jj+1,jk) - zprn(ji-1,jj+1,jk) ) * umask(ji-1,jj+1,jk) / e1u(ji-1,jj+1)   &
-                    + ( zprn(ji+1,jj+1,jk) - zprn(ji  ,jj+1,jk) ) * umask(ji  ,jj+1,jk) / e1u(ji  ,jj+1)   &
-                    + ( zprn(ji  ,jj  ,jk) - zprn(ji-1,jj  ,jk) ) * umask(ji-1,jj  ,jk) / e1u(ji-1,jj  )   &
-                    + ( zprn(ji+1,jj  ,jk) - zprn(ji  ,jj  ,jk) ) * umask(ji  ,jj  ,jk) / e1u(ji  ,jj  )
-               zphv = 1. / rau0 * zphv * zmsv * vmask(ji,jj,jk)
-
-               zmsu = 1. / MAX(  vmask(ji+1,jj  ,jk) + vmask(ji  ,jj  ,jk)   &
-                               + vmask(ji+1,jj-1,jk) + vmask(ji  ,jj-1,jk) , 1.  )
-               zphu = ( zprn(ji+1,jj+1,jk) - zprn(ji+1,jj  ,jk) ) * vmask(ji+1,jj  ,jk) / e2v(ji+1,jj  )   &
-                    + ( zprn(ji  ,jj+1,jk) - zprn(ji  ,jj  ,jk) ) * vmask(ji  ,jj  ,jk) / e2v(ji  ,jj  )   &
-                    + ( zprn(ji+1,jj  ,jk) - zprn(ji+1,jj-1,jk) ) * vmask(ji+1,jj-1,jk) / e2v(ji+1,jj-1)   &
-                    + ( zprn(ji  ,jj  ,jk) - zprn(ji  ,jj-1,jk) ) * vmask(ji  ,jj-1,jk) / e2v(ji  ,jj-1)
-               zphu = 1. / rau0 * zphu * zmsu * umask(ji,jj,jk)
-
-               ! Compute the geostrophic velocities
-               un(ji,jj,jk) = -2. * zphu / ( ff(ji,jj) + ff(ji  ,jj-1) )
-               vn(ji,jj,jk) =  2. * zphv / ( ff(ji,jj) + ff(ji-1,jj  ) )
-            END DO
-         END DO
-      END DO
-
-      IF(lwp) WRITE(numout,*) '         we force to zero bottom velocity'
-
-      ! Susbtract the bottom velocity (level jpk-1 for flat bottom case)
-      ! to have a zero bottom velocity
-
-      DO jk = 1, jpkm1
-         un(:,:,jk) = ( un(:,:,jk) - un(:,:,jpkm1) ) * umask(:,:,jk)
-         vn(:,:,jk) = ( vn(:,:,jk) - vn(:,:,jpkm1) ) * vmask(:,:,jk)
-      END DO
-
-      CALL lbc_lnk( un, 'U', -1. )
-      CALL lbc_lnk( vn, 'V', -1. )
-      
-      ub(:,:,:) = un(:,:,:)
-      vb(:,:,:) = vn(:,:,:)
-      
-      !
-!!gm  Check  here call to div_hor should not be necessary
-!!gm         div_hor call runoffs  not sure they are defined at that level
-      CALL div_hor( nit000 )            ! now horizontal divergence
-      !
-      CALL wrk_dealloc( jpi,jpj,jpk,   zprn)
-      !
-   END SUBROUTINE istate_uvg
-
-   !!=====================================================================
 END MODULE istate
Index: branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/usrdef.F90
===================================================================
--- branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/usrdef.F90	(revision 6484)
+++ branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/usrdef.F90	(revision 6579)
@@ -2,27 +2,17 @@
    !!==============================================================================
    !!                       ***  MODULE usrdef   ***
-   !! User defined module:  set analytically the minimum information needed to
-   !!       define a configuration (domain, initial state, forcing)
-   !!
-   !!              ===  GYRE configuration given as an example  ===
-   !! 
-   !!        This module is given as an example, it can be modified
-   !!      following the user's desiderata. 
-   !!      second order accuracy schemes.
+   !! User defined module: used like example to define init, sbc, bathy, ...
    !!==============================================================================
-   !! History :  NEMO ! 2016-03  (S. Flavoni)  Original code 
+   !! History :  NEMO ! 2016-03  (S. Flavoni) 
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
-   !!   usr_def_hgr       : compute the horizontal mesh 
-   !!   usr_def_zgr       : compute the vertical mesh 
-   !!  to be defined:
+   !!   usr_def_hgr       : initialize the horizontal mesh 
    !!   usr_def_ini       : initial state 
-   !!   usr_def_sbc       : compute the surface bounday conditions
-   !!   usr_def_xxx       : initialize the xxx
+   !!   usr_def_sbc       : initialize the surface bounday conditions
    !!----------------------------------------------------------------------
    USE step_oce       ! module used in the ocean time stepping module (step.F90)
-   USE mppini         ! shared/distributed memory setting (mpp_init routine)
-   USE phycst         ! physical constants
+   USE mppini         ! shared/distributed memory setting (mpp_init routine)             
+   USE phycst         ! physical constants 
    IMPLICIT NONE
    PRIVATE
@@ -30,20 +20,19 @@
 
    PUBLIC usr_def_hgr
-   PUBLIC usr_def_zgr
+   PUBLIC usr_def_ini
 
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.7 , NEMO Consortium (2016)
-   !! $Id: $ 
+   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
+   !! $Id:$ 
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 CONTAINS
 
-
-   SUBROUTINE usr_def_hgr( nbench, jp_cfg, kff, pff,       & ! Coriolis parameter  (if domain not on the sphere)
+   SUBROUTINE usr_def_hgr( nbench, jp_cfg, kff    , pff   , & ! Coriolis parameter  (if domain not on the sphere)
            &               pglamt, pglamu, pglamv , pglamf, & ! gridpoints position (required)
            &               pgphit, pgphiu, pgphiv , pgphif, & !         
            &               pe1t  , pe1u  , pe1v   , pe1f  , & ! scale factors (required)
            &               pe2t  , pe2u  , pe2v   , pe2f  , & !
-           &               ke1e2u_v                       )! u- & v-surfaces (if gridsize reduction is used in strait(s))
+           &               ke1e2u_v                       )   ! u- & v-surfaces (if gridsize reduction is used in strait(s))
       !!-------------------------------------------------------------------------------------------------
       !!                  ***  ROUTINE usr_def_hgr  ***
@@ -76,6 +65,6 @@
       !!-------------------------------------------------------------------------------
       INTEGER  ::   ji, jj               ! dummy loop indices
-      REAL(wp) ::   zlam1, zlam0, zcos_alpha, zim1 , zjm1 , ze1  , ze1deg, zf0   ! local scalars
-      REAL(wp) ::   zphi1, zphi0, zsin_alpha, zim05, zjm05, zbeta, znorme        ! local scalars
+      REAL(wp) ::   zlam1, zlam0, zcos_alpha, zim1 , zjm1 , ze1  , ze1deg, zf0 ! local scalars
+      REAL(wp) ::   zphi1, zphi0, zsin_alpha, zim05, zjm05, zbeta, znorme      ! local scalars
       !!-------------------------------------------------------------------------------
       !
@@ -91,80 +80,58 @@
       zphi1 =  29._wp
       ! resolution in meters
-      ze1 = 106000. / REAL( jp_cfg , wp )            
+      ze1 = 106000. / REAL( jp_cfg , wp )
       ! benchmark: forced the resolution to be about 100 km
-      IF( nbench /= 0 )   ze1 = 106000._wp   
-      zsin_alpha = - SQRT( 2._wp ) * 0.5_wp
-      zcos_alpha =   SQRT( 2._wp ) * 0.5_wp
-      ze1deg = ze1 / (ra * rad)
-      IF( nbench /= 0 )   ze1deg = ze1deg / REAL( jp_cfg , wp )   ! benchmark: keep the lat/+lon
-      !                                                           ! at the right jp_cfg resolution
-      zlam0 = zlam1 + zcos_alpha * ze1deg * REAL( jpjglo-2 , wp )
-      zphi0 = zphi1 + zsin_alpha * ze1deg * REAL( jpjglo-2 , wp )
-      !
-      IF( nprint==1 .AND. lwp )   THEN
-         WRITE(numout,*) '          ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha
-         WRITE(numout,*) '          ze1deg', ze1deg, 'zlam0', zlam0, 'zphi0', zphi0
-      ENDIF
-      !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zim1 = REAL( ji + nimpp - 1 ) - 1.   ;   zim05 = REAL( ji + nimpp - 1 ) - 1.5
-            zjm1 = REAL( jj + njmpp - 1 ) - 1.   ;   zjm05 = REAL( jj + njmpp - 1 ) - 1.5
-            !
-            !glamt(i,j) longitude at T-point
-            !gphit(i,j) latitude at T-point  
-            pglamt(ji,jj) = zlam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
-            pgphit(ji,jj) = zphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
-            !
-            !glamu(i,j) longitude at U-point
-            !gphiu(i,j) latitude at U-point
-            pglamu(ji,jj) = zlam0 + zim1  * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
-            pgphiu(ji,jj) = zphi0 - zim1  * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
-            !
-            !glamv(i,j) longitude at V-point
-            !gphiv(i,j) latitude at V-point
-            pglamv(ji,jj) = zlam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
-            pgphiv(ji,jj) = zphi0 - zim05 * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
-            !
-            !glamf(i,j) longitude at F-point
-            !gphif(i,j) latitude at F-point 
-            pglamf(ji,jj) = zlam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
-            pgphif(ji,jj) = zphi0 - zim1  * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
-         END DO
-      END DO
-      !
-      !       !== Horizontal scale factors ==! (in meters)
-      !                     
-      !                             ! constant grid spacing
-      pe1t(:,:) =  ze1     ;      pe2t(:,:) = ze1
-      pe1u(:,:) =  ze1     ;      pe2u(:,:) = ze1
-      pe1v(:,:) =  ze1     ;      pe2v(:,:) = ze1
-      pe1f(:,:) =  ze1     ;      pe2f(:,:) = ze1
-      !                             ! NO reduction of grid size in some straits 
-      ke1e2u_v = 0                  !    ==>> u_ & v_surfaces will be computed in dom_ghr routine
-      !
-      !
-
-      !                      !==  Coriolis parameter  ==!
-      kff = 1                                                           !  indicate not to compute ff afterward
-      !
-      zbeta = 2. * omega * COS( rad * zphi1 ) / ra                      ! beta at latitude zphi1
-      !SF we overwrite zphi0 (south point in latitude) used just above to define pgphif (value of zphi0=15.5190567531966)
-      !SF for computation of coriolis we keep the parameter of XXXX
-      zphi0 = 15._wp                                                     !  latitude of the most southern grid point  
-      zf0   = 2. * omega * SIN( rad * zphi0 )                            !  compute f0 1st point south
-      !
-      pff(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south)
-      !
-      IF(lwp) WRITE(numout,*) '                           beta-plane used. beta = ', zbeta, ' 1/(s.m)'
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('usr_def_hgr')
       !
    END SUBROUTINE usr_def_hgr
-  
+
+
    SUBROUTINE usr_def_zgr()           ! Coriolis parameter  (if domain not on the sphere)
        ! subroutine for vertical grid
    END SUBROUTINE usr_def_zgr
 
+ 
+   SUBROUTINE usr_def_ini( pts )
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE usr_def_ini  ***
+      !! 
+      !! ** Purpose :   Initialization of the dynamics and tracers
+      !!                Here GYRE configuration example : (double gyre with rotated domain)
+      !!
+      !! ** Method  : - set temprature field
+      !!              - set salinity   field
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   ptsd   ! T & S data
+      !
+      INTEGER :: ji, jj, jk  ! dummy loop indices
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'usr_def_ini : analytical definition of initial state '
+      IF(lwp) WRITE(numout,*) '              T and S profiles deduced from LEVITUS '
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+      !
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               pts(ji,jj,jk,jp_tem) =  (  (  16. - 12. * TANH( (gdept_n(ji,jj,jk) - 400) / 700 ) )   &
+                    &           * (-TANH( (500-gdept_n(ji,jj,jk)) / 150 ) + 1) / 2               &
+                    &           + ( 15. * ( 1. - TANH( (gdept_n(ji,jj,jk)-50.) / 1500.) )        &
+                    &           - 1.4 * TANH((gdept_n(ji,jj,jk)-100.) / 100.)                    &
+                    &           + 7.  * (1500. - gdept_n(ji,jj,jk)) / 1500.)                     &
+                    &           * (-TANH( (gdept_n(ji,jj,jk) - 500) / 150) + 1) / 2  ) * tmask(ji,jj,jk)
+
+               pts(ji,jj,jk,jp_sal) =  (  (  36.25 - 1.13 * TANH( (gdept_n(ji,jj,jk) - 305) / 460 ) )  &
+                    &         * (-TANH((500 - gdept_n(ji,jj,jk)) / 150) + 1) / 2              &
+                    &         + ( 35.55 + 1.25 * (5000. - gdept_n(ji,jj,jk)) / 5000.          &
+                    &         - 1.62 * TANH( (gdept_n(ji,jj,jk) - 60.  ) / 650. )             &
+                    &         + 0.2  * TANH( (gdept_n(ji,jj,jk) - 35.  ) / 100. )             &
+                    &         + 0.2  * TANH( (gdept_n(ji,jj,jk) - 1000.) / 5000.) )           &
+                    &         * (-TANH((gdept_n(ji,jj,jk) - 500) / 150) + 1) / 2  ) * tmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !   
+   END SUBROUTINE usr_def_ini
+
    !!======================================================================
 END MODULE usrdef
