Index: /branches/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src/create_meshmask.f90
===================================================================
--- /branches/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src/create_meshmask.f90	(revision 7232)
+++ /branches/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src/create_meshmask.f90	(revision 7233)
@@ -66,9 +66,11 @@
 !>
 !>    * _input files namelist (namin)_:<br/>
-!>       - cn_bathy  : Bathymetry file
-!>       - cn_coord  : coordinate file (in_mshhgr=0)
-!>       - cn_isfdep : Iceshelf draft  (ln_isfcav=true)
-!>       - in_perio  : NEMO periodicity
-!>       - ln_closea :
+!>       - cn_bathy     : Bathymetry file
+!>       - cn_varbathy  : Bathymetry variable name
+!>       - cn_coord     : coordinate file (in_mshhgr=0)
+!>       - cn_isfdep    : Iceshelf draft  (ln_isfcav=true)
+!>       - cn_varisfdep : Iceshelf draft variable name (ln_isfcav=true)
+!>       - in_perio     : NEMO periodicity
+!>       - ln_closea    :
 !>
 !>    * _horizontal grid namelist (namhgr)_:<br/>
@@ -164,4 +166,5 @@
 !       - ln_zoom   : use zoom (namzoom)
 !>       - ln_c1d    : use configuration 1D
+!>       - ln_e3_dep : vertical scale factors =T: e3.=dk[depth] =F: old definition
 !
 !    * _zoom namelist (namzoom)_:<br/>
@@ -175,4 +178,5 @@
 !
 !>    * _output namelist (namout)_:<br/>
+!>       - cn_domcfg    : output file name
 !>       - in_msh       : number of output file and contain (0-9)
 !>       - in_nproc     : number of processor to be used
@@ -204,4 +208,6 @@
 !> - do not use anymore special case for ORCA grid
 !> - allow to write domain_cfg file
+!> @date November, 2016
+!> - choose vertical scale factors (e3.=dk[depth] or old definition)
 !>
 !> @note Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -298,6 +304,8 @@
    ! namin
    CHARACTER(LEN=lc) :: cn_bathy    = ''
+   CHARACTER(LEN=lc) :: cn_varbathy = ''
    CHARACTER(LEN=lc) :: cn_coord    = ''
    CHARACTER(LEN=lc) :: cn_isfdep   = ''
+   CHARACTER(LEN=lc) :: cn_varisfdep= 'isfdraft'
    INTEGER(i4)       :: in_perio    = -1
    LOGICAL           :: ln_closea   = .TRUE.
@@ -313,7 +321,8 @@
 
    ! namlbc
-   REAL(sp)          :: rn_shlat    = 2.
+   REAL(sp)          :: rn_shlat    = -100. !2.
 
    ! namout
+   CHARACTER(LEN=lc) :: cn_domcfg   = 'domain_cfg.nc'
    INTEGER(i4)       :: in_msh      = 0 
    CHARACTER(LEN=lc) :: cn_type     = 'cdf'
@@ -334,6 +343,8 @@
    NAMELIST /namin/  &  !< input namelist
    &  cn_bathy,      &  !< Bathymetry file
+   &  cn_varbathy,   &  !< Bathymetry variable name
    &  cn_coord,      &  !< Coordinate file (in_mshhgr=0)
    &  cn_isfdep,     &  !< Iceshelf draft  (ln_isfcav=true)
+   &  cn_varisfdep,  &  !< Iceshelf draft variable name (ln_isfcav=true)
    &  in_perio,      &  !< NEMO periodicity
    &  ln_closea
@@ -352,4 +363,5 @@
 
    NAMELIST /namout/ &  !< output namlist
+   &  cn_domcfg,     &  !< output file name
    &  in_msh,        &  !< number of output file (1,2,3)
    &  in_nproc,      &  !< number of processor to be used
@@ -449,16 +461,18 @@
    WRITE(*,*) 'DIMENSION TO BE USED :',jpi,jpj,jpk
 
-   IF( ln_zco )   THEN
-      WRITE(*,*) 'VARIABLE READ : Bathy_level'
-      tl_bathy=iom_mpp_read_var(tl_mpp, 'Bathy_level')
-   ELSEIF( ln_zps .OR. ln_sco )   THEN
-      IF ( ln_isfcav ) THEN
-         WRITE(*,*) 'VARIABLE READ : Bathymetry_isf'
-         tl_bathy=iom_mpp_read_var(tl_mpp, 'Bathymetry_isf' )
-      ELSE
-         WRITE(*,*) 'VARIABLE READ : Bathymetry'
-         tl_bathy=iom_mpp_read_var(tl_mpp, 'Bathymetry' )
-      END IF
+   ! read variable
+   IF( TRIM(cn_varbathy) == '' )THEN
+      IF( ln_zco )THEN
+         cn_varbathy='Bathy_level'
+      ELSEIF( ln_zps .OR. ln_sco )THEN
+         IF( ln_isfcav )THEN
+            cn_varbathy='Bathymetry_isf'
+         ELSE
+            cn_varbathy='Bathymetry'
+         ENDIF
+      ENDIF
    ENDIF
+   WRITE(*,*) 'VARIABLE READ : '//TRIM(cn_varbathy)
+   tl_bathy=iom_mpp_read_var(tl_mpp, TRIM(cn_varbathy))
    CALL iom_mpp_close(tl_mpp)
 
@@ -469,4 +483,6 @@
  
    IF ( ln_isfcav ) THEN
+      WRITE(*,*) 'ICESHELF DRAFT FILE TO BE USED:',TRIM(cn_isfdep)
+      WRITE(*,*) 'ICESHELF VARIABLE READ : '//TRIM(cn_varisfdep)
       ! open Iceshelf draft
       IF( cn_isfdep /= '' )THEN
@@ -481,5 +497,5 @@
       CALL iom_mpp_open(tl_mpp)
       IF( ln_zps .OR. ln_sco )   THEN
-         tl_risfdep=iom_mpp_read_var(tl_mpp, 'isf_draft')
+         tl_risfdep=iom_mpp_read_var(tl_mpp, cn_varisfdep)
       ENDIF
       CALL iom_mpp_close(tl_mpp)
@@ -488,5 +504,5 @@
       dl_tmp2D(:,:)=0._dp
 
-      tl_risfdep=var_init('isf_draft',dl_tmp2D(:,:), id_type=NF90_FLOAT)
+      tl_risfdep=var_init(cn_varisfdep, dl_tmp2D(:,:), id_type=NF90_DOUBLE)
 
       DEALLOCATE(dl_tmp2D)
@@ -561,5 +577,5 @@
          !                                  ! create 'domain_cfg.nc' file
          !                                  ! ============================
-         tl_mppout0=mpp_init( 'domain_cfg', tg_tmask, &
+         tl_mppout0=mpp_init( cn_domcfg, tg_tmask, &
          &                    in_niproc, in_njproc, in_nproc, &
          &                    cd_type=cn_type )
@@ -789,5 +805,7 @@
    &                            tg_mbkt%d_value(:,:,:,:)
    !
-   IF( ll_domcfg ) tg_mbathy%c_name='bottom_level'
+   IF( ll_domcfg )THEN
+      tg_mbathy%c_name='bottom_level'
+   ENDIF
    CALL mpp_add_var(tl_mppzgr, tg_mbathy)
    CALL var_clean(tg_mbathy)
@@ -862,7 +880,16 @@
       CALL mpp_add_var(tl_mppzgr, tg_e3v_0)
       CALL var_clean(tg_e3v_0)
+      ! e3f_0
+      CALL mpp_add_var(tl_mppzgr, tg_e3f_0)
+      CALL var_clean(tg_e3f_0)
       ! e3w_0
       CALL mpp_add_var(tl_mppzgr, tg_e3w_0)
       CALL var_clean(tg_e3w_0)
+      ! e3uw_0
+      CALL mpp_add_var(tl_mppzgr, tg_e3uw_0)
+      CALL var_clean(tg_e3uw_0)
+      ! e3vw_0
+      CALL mpp_add_var(tl_mppzgr, tg_e3vw_0)
+      CALL var_clean(tg_e3vw_0)
 
       ! Max. grid stiffness ratio
@@ -873,17 +900,19 @@
 
       ! stretched system
-      ! gdept_1d
-      CALL mpp_add_var(tl_mppzgr, tg_gdept_1d)
-      CALL var_clean(tg_gdept_1d)
-      ! gdepw_1d
-      CALL mpp_add_var(tl_mppzgr, tg_gdepw_1d)
-      CALL var_clean(tg_gdepw_1d)
-      
-      ! gdept_0
-      CALL mpp_add_var(tl_mppzgr, tg_gdept_0)      
-      CALL var_clean(tg_gdept_0)
-      ! gdepw_0
-      CALL mpp_add_var(tl_mppzgr, tg_gdepw_0)      
-      CALL var_clean(tg_gdepw_0)
+      IF( .NOT. tl_namz%l_e3_dep )THEN
+         ! gdept_1d
+         CALL mpp_add_var(tl_mppzgr, tg_gdept_1d)
+         CALL var_clean(tg_gdept_1d)
+         ! gdepw_1d
+         CALL mpp_add_var(tl_mppzgr, tg_gdepw_1d)
+         CALL var_clean(tg_gdepw_1d)
+         
+         ! gdept_0
+         CALL mpp_add_var(tl_mppzgr, tg_gdept_0)      
+         CALL var_clean(tg_gdept_0)
+         ! gdepw_0
+         CALL mpp_add_var(tl_mppzgr, tg_gdepw_0)      
+         CALL var_clean(tg_gdepw_0)
+      ENDIF
 
    ENDIF
@@ -928,44 +957,45 @@
       IF( ll_domcfg .OR. in_msh <= 3 ) THEN ! 3D depth
          
-         ! gdept_0
-         CALL mpp_add_var(tl_mppzgr, tg_gdept_0)
+         IF( .NOT. tl_namz%l_e3_dep )THEN
       
-         ! gdepu, gdepv
-         IF( .NOT. ll_domcfg )THEN
-            ALLOCATE(dl_tmp3D(jpi,jpj,jpk))
-            dl_tmp3D(:,:,:)=dp_fill
-
-            tl_gdepu=var_init('gdepu',dl_tmp3D(:,:,:), id_type=NF90_FLOAT)
-            tl_gdepv=var_init('gdepv',dl_tmp3D(:,:,:), id_type=NF90_FLOAT)
-
-            DEALLOCATE(dl_tmp3D)
-            DO jk = 1,jpk   
-               DO jj = 1, jpj-1   
-                  DO ji = 1, jpi-1   ! vector opt.
-                     tl_gdepu%d_value(ji,jj,jk,1) = MIN( tg_gdept_0%d_value(ji  ,jj  ,jk,1) , &
-                        &                                tg_gdept_0%d_value(ji+1,jj  ,jk,1) )
-
-                     tl_gdepv%d_value(ji,jj,jk,1) = MIN( tg_gdept_0%d_value(ji  ,jj  ,jk,1) , &
-                        &                                tg_gdept_0%d_value(ji  ,jj+1,jk,1) )
+            ! gdepu, gdepv
+            IF( .NOT. ll_domcfg )THEN
+               ALLOCATE(dl_tmp3D(jpi,jpj,jpk))
+               dl_tmp3D(:,:,:)=dp_fill
+
+               tl_gdepu=var_init('gdepu',dl_tmp3D(:,:,:), id_type=NF90_FLOAT)
+               tl_gdepv=var_init('gdepv',dl_tmp3D(:,:,:), id_type=NF90_FLOAT)
+
+               DEALLOCATE(dl_tmp3D)
+               DO jk = 1,jpk   
+                  DO jj = 1, jpj-1   
+                     DO ji = 1, jpi-1   ! vector opt.
+                        tl_gdepu%d_value(ji,jj,jk,1) = MIN( tg_gdept_0%d_value(ji  ,jj  ,jk,1) , &
+                           &                                tg_gdept_0%d_value(ji+1,jj  ,jk,1) )
+
+                        tl_gdepv%d_value(ji,jj,jk,1) = MIN( tg_gdept_0%d_value(ji  ,jj  ,jk,1) , &
+                           &                                tg_gdept_0%d_value(ji  ,jj+1,jk,1) )
+                     END DO   
                   END DO   
-               END DO   
-            END DO         
-            CALL lbc_lnk( tl_gdepu%d_value(:,:,:,1), 'U', in_perio, 1._dp )
-            CALL lbc_lnk( tl_gdepv%d_value(:,:,:,1), 'V', in_perio, 1._dp )
-
-            ! gdepu
-            CALL mpp_add_var(tl_mppzgr, tl_gdepu)
-            CALL var_clean(tl_gdepu)
-            ! gdepv
-            CALL mpp_add_var(tl_mppzgr, tl_gdepv)
-            CALL var_clean(tl_gdepv)
+               END DO         
+               CALL lbc_lnk( tl_gdepu%d_value(:,:,:,1), 'U', in_perio, 1._dp )
+               CALL lbc_lnk( tl_gdepv%d_value(:,:,:,1), 'V', in_perio, 1._dp )
+
+               ! gdepu
+               CALL mpp_add_var(tl_mppzgr, tl_gdepu)
+               CALL var_clean(tl_gdepu)
+               ! gdepv
+               CALL mpp_add_var(tl_mppzgr, tl_gdepv)
+               CALL var_clean(tl_gdepv)
+            ENDIF
+
+            ! gdept_0
+            CALL mpp_add_var(tl_mppzgr, tg_gdept_0)
+            CALL var_clean(tg_gdept_0)
+
+            ! gdepw_0
+            CALL mpp_add_var(tl_mppzgr, tg_gdepw_0)
+            CALL var_clean(tg_gdepw_0)
          ENDIF
-
-         ! clean
-         CALL var_clean(tg_gdept_0)
-
-         ! gdepw_0
-         CALL mpp_add_var(tl_mppzgr, tg_gdepw_0)
-         CALL var_clean(tg_gdepw_0)
 
       ELSE ! 2D bottom depth
@@ -998,13 +1028,6 @@
    ENDIF
 
-   IF( ln_zps .OR. ln_zco )THEN ! z-coordinate 
-      ! depth
-      ! gdept_1d
-      CALL mpp_add_var(tl_mppzgr, tg_gdept_1d)
-      CALL var_clean(tg_gdept_1d)
-      ! gdepw_1d
-      CALL mpp_add_var(tl_mppzgr, tg_gdepw_1d)
-      CALL var_clean(tg_gdepw_1d)
-      ! scale factors
+   ! scale factors
+   IF( ll_domcfg )THEN
       ! e3t_1d
       CALL mpp_add_var(tl_mppzgr, tg_e3t_1d)
@@ -1015,4 +1038,16 @@
    ENDIF
 
+   IF( ln_zps .OR. ln_zco )THEN ! z-coordinate 
+      IF( .NOT. tl_namz%l_e3_dep )THEN
+         ! depth
+         ! gdept_1d
+         CALL mpp_add_var(tl_mppzgr, tg_gdept_1d)
+         CALL var_clean(tg_gdept_1d)
+         ! gdepw_1d
+         CALL mpp_add_var(tl_mppzgr, tg_gdepw_1d)
+         CALL var_clean(tg_gdepw_1d)
+      ENDIF
+   ENDIF
+
    ! define global attributes
    ALLOCATE(tl_gatt(ip_maxatt))
@@ -1021,4 +1056,5 @@
 
 
+   IF( in_msh == 0 ) in_msh=1
    SELECT CASE ( MOD(in_msh, 3) )
       CASE ( 1 )
@@ -1404,5 +1440,4 @@
 !      ENDIF
 
-      IF( .NOT. ld_domcfg )THEN
       ! 2. Ocean/land mask at u-,  v-, and z-points (computed from tmask)
       ! -------------------------------------------
@@ -1417,8 +1452,10 @@
             !END DO
             !DO ji = 1, jpi-1  ! NO vector opt.
-               tg_fmask%d_value(ji,jj,jk,1) = tg_tmask%d_value(ji  ,jj  ,jk,1) * &
-                                            & tg_tmask%d_value(ji+1,jj  ,jk,1) * &
-                                            & tg_tmask%d_value(ji  ,jj+1,jk,1) * &
-                                            & tg_tmask%d_value(ji+1,jj+1,jk,1)
+               IF( .NOT. ld_domcfg )THEN
+                  tg_fmask%d_value(ji,jj,jk,1) = tg_tmask%d_value(ji  ,jj  ,jk,1) * &
+                                               & tg_tmask%d_value(ji+1,jj  ,jk,1) * &
+                                               & tg_tmask%d_value(ji  ,jj+1,jk,1) * &
+                                               & tg_tmask%d_value(ji+1,jj+1,jk,1)
+               ENDIF
             ENDDO
          ENDDO
@@ -1445,5 +1482,7 @@
       CALL lbc_lnk( tg_umask%d_value  (:,:,:,1), 'U', td_nam%i_perio, 1._dp )      ! Lateral boundary conditions
       CALL lbc_lnk( tg_vmask%d_value  (:,:,:,1), 'V', td_nam%i_perio, 1._dp )
-      CALL lbc_lnk( tg_fmask%d_value  (:,:,:,1), 'F', td_nam%i_perio, 1._dp )
+      IF( .NOT. ld_domcfg )THEN
+         CALL lbc_lnk( tg_fmask%d_value  (:,:,:,1), 'F', td_nam%i_perio, 1._dp )
+      ENDIF
 !      CALL lbc_lnk( tg_ssumask%d_value(:,:,:,1), 'U', td_nam%i_perio, 1._dp )      ! Lateral boundary conditions
 !      CALL lbc_lnk( tg_ssvmask%d_value(:,:,:,1), 'V', td_nam%i_perio, 1._dp )
@@ -1467,38 +1506,39 @@
       ! Lateral boundary conditions on velocity (modify fmask)
       ! ---------------------------------------     
-      ALLOCATE( zwf(jpi,jpj) )
-      DO jk = 1, jpk
-         zwf(:,:) = tg_fmask%d_value(:,:,jk,1)         
-         DO jj = 2, jpj-1
-            DO ji = 2, jpi-1   ! vector opt.
-               IF( tg_fmask%d_value(ji,jj,jk,1) == 0._dp )THEN
-                  tg_fmask%d_value(ji,jj,jk,1) = rn_shlat * &
-                                               & MIN(1._dp , MAX(zwf(ji+1,jj), zwf(ji,jj+1), &
-                                               &                 zwf(ji-1,jj), zwf(ji,jj-1)) )
+      IF( .NOT. ld_domcfg )THEN
+         ALLOCATE( zwf(jpi,jpj) )
+         DO jk = 1, jpk
+            zwf(:,:) = tg_fmask%d_value(:,:,jk,1)         
+            DO jj = 2, jpj-1
+               DO ji = 2, jpi-1   ! vector opt.
+                  IF( tg_fmask%d_value(ji,jj,jk,1) == 0._dp )THEN
+                     tg_fmask%d_value(ji,jj,jk,1) = rn_shlat * &
+                                                  & MIN(1._dp , MAX(zwf(ji+1,jj), zwf(ji,jj+1), &
+                                                  &                 zwf(ji-1,jj), zwf(ji,jj-1)) )
+                  ENDIF
+               END DO
+            END DO
+            DO jj = 2, jpj-1
+               IF( tg_fmask%d_value(1,jj,jk,1) == 0._dp )THEN
+                  tg_fmask%d_value(1  ,jj,jk,1) = rn_shlat * &
+                                                & MIN(1._dp, MAX(zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1)))
+               ENDIF
+               IF( tg_fmask%d_value(jpi,jj,jk,1) == 0._dp )THEN
+                  tg_fmask%d_value(jpi,jj,jk,1) = rn_shlat * &
+                                                & MIN(1._wp, MAX(zwf(jpi,jj+1), zwf(jpi-1,jj), zwf(jpi,jj-1)))
+               ENDIF
+            END DO         
+            DO ji = 2, jpi-1
+               IF( tg_fmask%d_value(ji,1,jk,1) == 0._dp )THEN
+                  tg_fmask%d_value(ji, 1 ,jk,1) = rn_shlat * &
+                                                & MIN(1._dp, MAX(zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1)))
+               ENDIF
+               IF( tg_fmask%d_value(ji,jpj,jk,1) == 0._dp )THEN
+                  tg_fmask%d_value(ji,jpj,jk,1) = rn_shlat * &
+                                                & MIN(1._dp, MAX(zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpj-1)))
                ENDIF
             END DO
          END DO
-         DO jj = 2, jpj-1
-            IF( tg_fmask%d_value(1,jj,jk,1) == 0._dp )THEN
-               tg_fmask%d_value(1  ,jj,jk,1) = rn_shlat * &
-                                             & MIN(1._dp, MAX(zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1)))
-            ENDIF
-            IF( tg_fmask%d_value(jpi,jj,jk,1) == 0._dp )THEN
-               tg_fmask%d_value(jpi,jj,jk,1) = rn_shlat * &
-                                             & MIN(1._wp, MAX(zwf(jpi,jj+1), zwf(jpi-1,jj), zwf(jpi,jj-1)))
-            ENDIF
-         END DO         
-         DO ji = 2, jpi-1
-            IF( tg_fmask%d_value(ji,1,jk,1) == 0._dp )THEN
-               tg_fmask%d_value(ji, 1 ,jk,1) = rn_shlat * &
-                                             & MIN(1._dp, MAX(zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1)))
-            ENDIF
-            IF( tg_fmask%d_value(ji,jpj,jk,1) == 0._dp )THEN
-               tg_fmask%d_value(ji,jpj,jk,1) = rn_shlat * &
-                                             & MIN(1._dp, MAX(zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpj-1)))
-            ENDIF
-         END DO
-      END DO
-      DEALLOCATE( zwf )
+         DEALLOCATE( zwf )
 
 !      IF( td_nam%c_cfg == "orca" .AND. td_nam%i_cfg == 2 )THEN   ! ORCA_R2 configuration
@@ -1576,8 +1616,8 @@
 !      ENDIF
       !
-      CALL lbc_lnk( tg_fmask%d_value(:,:,:,1), 'F', td_nam%i_perio, 1._dp )      ! Lateral boundary conditions on fmask
+         CALL lbc_lnk( tg_fmask%d_value(:,:,:,1), 'F', td_nam%i_perio, 1._dp )      ! Lateral boundary conditions on fmask
+      ENDIF
 
 !      DEALLOCATE( imsk )
-      ENDIF ! ld_domcfg
 
 !      DEALLOCATE( dl_tpol )
@@ -1700,5 +1740,5 @@
             ji=ji+1 ; td_att(ji)=att_init("bb",td_namz%d_bb)
          ELSEIF( td_namz%l_s_sf12 )THEN
-            IF( td_namz%l_sigcrit ) ji=ji+1 ; td_att(ji)=att_init("sigma below critical depth","activated")
+            IF( td_namz%l_sigcrit ) ji=ji+1 ; td_att(ji)=att_init("sigma_below_critical_depth","activated")
             ji=ji+1 ; td_att(ji)=att_init("alpha",td_namz%d_alpha)
             ji=ji+1 ; td_att(ji)=att_init("efold",td_namz%d_efold)
Index: /branches/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src/grid_hgr.f90
===================================================================
--- /branches/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src/grid_hgr.f90	(revision 7232)
+++ /branches/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src/grid_hgr.f90	(revision 7233)
@@ -273,17 +273,15 @@
 !      tg_ssfmask  = var_init('ssfmask',dl_tmp2D(:,:)  , dd_fill=dp_fill_i1, id_type=NF90_BYTE)
 
-      dl_tmp2D(:,:)=dp_fill_sp
-
-      tg_glamt    = var_init('glamt',dl_tmp2D(:,:)   , dd_fill=dp_fill_sp, id_type=NF90_FLOAT)
-      tg_glamu    = var_init('glamu',dl_tmp2D(:,:)   , dd_fill=dp_fill_sp, id_type=NF90_FLOAT)
-      tg_glamv    = var_init('glamv',dl_tmp2D(:,:)   , dd_fill=dp_fill_sp, id_type=NF90_FLOAT)
-      tg_glamf    = var_init('glamf',dl_tmp2D(:,:)   , dd_fill=dp_fill_sp, id_type=NF90_FLOAT)
-
-      tg_gphit    = var_init('gphit',dl_tmp2D(:,:)   , dd_fill=dp_fill_sp, id_type=NF90_FLOAT)
-      tg_gphiu    = var_init('gphiu',dl_tmp2D(:,:)   , dd_fill=dp_fill_sp, id_type=NF90_FLOAT)
-      tg_gphiv    = var_init('gphiv',dl_tmp2D(:,:)   , dd_fill=dp_fill_sp, id_type=NF90_FLOAT)
-      tg_gphif    = var_init('gphif',dl_tmp2D(:,:)   , dd_fill=dp_fill_sp, id_type=NF90_FLOAT)
-
       dl_tmp2D(:,:)=dp_fill
+
+      tg_glamt    = var_init('glamt',dl_tmp2D(:,:)   , dd_fill=dp_fill, id_type=NF90_DOUBLE)
+      tg_glamu    = var_init('glamu',dl_tmp2D(:,:)   , dd_fill=dp_fill, id_type=NF90_DOUBLE)
+      tg_glamv    = var_init('glamv',dl_tmp2D(:,:)   , dd_fill=dp_fill, id_type=NF90_DOUBLE)
+      tg_glamf    = var_init('glamf',dl_tmp2D(:,:)   , dd_fill=dp_fill, id_type=NF90_DOUBLE)
+
+      tg_gphit    = var_init('gphit',dl_tmp2D(:,:)   , dd_fill=dp_fill, id_type=NF90_DOUBLE)
+      tg_gphiu    = var_init('gphiu',dl_tmp2D(:,:)   , dd_fill=dp_fill, id_type=NF90_DOUBLE)
+      tg_gphiv    = var_init('gphiv',dl_tmp2D(:,:)   , dd_fill=dp_fill, id_type=NF90_DOUBLE)
+      tg_gphif    = var_init('gphif',dl_tmp2D(:,:)   , dd_fill=dp_fill, id_type=NF90_DOUBLE)
 
       tg_e1t      = var_init('e1t'  ,dl_tmp2D(:,:)   , dd_fill=dp_fill, id_type=NF90_DOUBLE)
@@ -316,7 +314,7 @@
 
       tg_tmask   = var_init('tmask'  ,dl_tmp3D(:,:,:), dd_fill=dp_fill_i1, id_type=NF90_BYTE)
+      tg_umask   = var_init('umask'  ,dl_tmp3D(:,:,:), dd_fill=dp_fill_i1, id_type=NF90_BYTE)
+      tg_vmask   = var_init('vmask'  ,dl_tmp3D(:,:,:), dd_fill=dp_fill_i1, id_type=NF90_BYTE)
       IF( .NOT. ld_domcfg )THEN
-         tg_umask   = var_init('umask'  ,dl_tmp3D(:,:,:), dd_fill=dp_fill_i1, id_type=NF90_BYTE)
-         tg_vmask   = var_init('vmask'  ,dl_tmp3D(:,:,:), dd_fill=dp_fill_i1, id_type=NF90_BYTE)
          tg_fmask   = var_init('fmask'  ,dl_tmp3D(:,:,:), dd_fill=dp_fill_i1, id_type=NF90_BYTE)
       ENDIF
@@ -380,7 +378,7 @@
 
       CALL var_clean(tg_tmask   )
+      CALL var_clean(tg_umask   )
+      CALL var_clean(tg_vmask   )
       IF( .NOT. ld_domcfg )THEN
-         CALL var_clean(tg_umask   )
-         CALL var_clean(tg_vmask   )
          CALL var_clean(tg_fmask   )
       ENDIF
@@ -533,5 +531,5 @@
       ! loop indices
       !----------------------------------------------------------------
-      CALL logger_info('GRID HGR FILL : define the horizontal mesh from ithe'//&
+      CALL logger_info('GRID HGR FILL : define the horizontal mesh from the'//&
       &  ' type of horizontal mesh mshhgr = '//TRIM(fct_str(td_nam%i_mshhgr)))
       IF( td_nam%i_mshhgr == 1 )THEN
@@ -657,13 +655,13 @@
 
       ! force output type
-      tg_glamt%i_type=NF90_FLOAT
-      tg_glamu%i_type=NF90_FLOAT
-      tg_glamv%i_type=NF90_FLOAT
-      tg_glamf%i_type=NF90_FLOAT
-
-      tg_gphit%i_type=NF90_FLOAT
-      tg_gphiu%i_type=NF90_FLOAT
-      tg_gphiv%i_type=NF90_FLOAT
-      tg_gphif%i_type=NF90_FLOAT
+      tg_glamt%i_type=NF90_DOUBLE
+      tg_glamu%i_type=NF90_DOUBLE
+      tg_glamv%i_type=NF90_DOUBLE
+      tg_glamf%i_type=NF90_DOUBLE
+
+      tg_gphit%i_type=NF90_DOUBLE
+      tg_gphiu%i_type=NF90_DOUBLE
+      tg_gphiv%i_type=NF90_DOUBLE
+      tg_gphif%i_type=NF90_DOUBLE
 
       tg_e1t  =iom_mpp_read_var(tl_coord, 'e1t')
Index: /branches/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src/grid_zgr.f90
===================================================================
--- /branches/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src/grid_zgr.f90	(revision 7232)
+++ /branches/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src/grid_zgr.f90	(revision 7233)
@@ -55,4 +55,6 @@
 !> - J, Paul : update from trunk (revision 6961): add wetting and drying, ice sheet coupling..
 !> - J, Paul : do not use anymore special case for ORCA grid.
+!> @date November, 2016
+!> - J, Paul : vertical scale factors e3. = dk[gdep] or old definition
 !>
 !> @note Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -99,7 +101,7 @@
    PUBLIC :: tg_e3v_0
    PUBLIC :: tg_e3w_0
-!   PUBLIC :: tg_e3f_0     !useless to create meshmask
-!   PUBLIC :: tg_e3uw_0    !useless to create meshmask
-!   PUBLIC :: tg_e3vw_0    !useless to create meshmask
+   PUBLIC :: tg_e3f_0     !useless to create meshmask
+   PUBLIC :: tg_e3uw_0    !useless to create meshmask
+   PUBLIC :: tg_e3vw_0    !useless to create meshmask
    
    PUBLIC :: tg_mbkt
@@ -144,5 +146,5 @@
    PRIVATE :: grid_zgr__isf_fill
 !   PRIVATE :: grid_zgr__isf_fill_e3x
-!   PRIVATE :: grid_zgr__isf_fill_e3uw
+   PRIVATE :: grid_zgr__isf_fill_e3uw
 !   PRIVATE :: grid_zgr__isf_fill_gdep3w_0
    PRIVATE :: grid_zgr__sco_fill
@@ -155,4 +157,5 @@
 
    TYPE TNAMZ
+
       CHARACTER(LEN=lc) :: c_coord   
       INTEGER(i4)       :: i_perio   
@@ -214,4 +217,5 @@
 !      LOGICAL           :: l_zoom
       LOGICAL           :: l_c1d
+      LOGICAL           :: l_e3_dep
                 
 !      CHARACTER(LEN=lc) :: c_cfz      
@@ -244,7 +248,7 @@
    TYPE(TVAR), SAVE :: tg_e3v_0     !      zps & sco
    TYPE(TVAR), SAVE :: tg_e3w_0     !      zps & sco
-   !TYPE(TVAR), SAVE :: tg_e3f_0
-   !TYPE(TVAR), SAVE :: tg_e3uw_0
-   !TYPE(TVAR), SAVE :: tg_e3vw_0
+   TYPE(TVAR), SAVE :: tg_e3f_0
+   TYPE(TVAR), SAVE :: tg_e3uw_0
+   TYPE(TVAR), SAVE :: tg_e3vw_0
    
    TYPE(TVAR), SAVE :: tg_mbkt      !zco & zps & sco
@@ -315,7 +319,4 @@
       dl_tmp2D(:,:)  =dp_fill_i2
 
-      tg_mbathy  =var_init('mbathy  ',dl_tmp2D(:,:)  , dd_fill=dp_fill_i2, id_type=NF90_SHORT)
-      tg_misfdep =var_init('misfdep ',dl_tmp2D(:,:)  , dd_fill=dp_fill_i2, id_type=NF90_SHORT)
-
       tg_mbkt    =var_init('mbkt    ',dl_tmp2D(:,:)  , dd_fill=dp_fill_i2, id_type=NF90_SHORT)
       !tg_mbku    =var_init('mbku    ',dl_tmp2D(:,:)  , dd_fill=dp_fill_i2, id_type=NF90_SHORT)
@@ -326,4 +327,9 @@
       !tg_mikf    =var_init('mikf    ',dl_tmp2D(:,:)  , dd_fill=dp_fill_i2, id_type=NF90_SHORT)
 
+      dl_tmp2D(:,:)  =dp_fill_i4
+
+      tg_mbathy  =var_init('mbathy  ',dl_tmp2D(:,:)  , dd_fill=dp_fill_i4, id_type=NF90_INT)
+      tg_misfdep =var_init('misfdep ',dl_tmp2D(:,:)  , dd_fill=dp_fill_i4, id_type=NF90_INT)
+
       dl_tmp2D(:,:)  =dp_fill
 
@@ -337,19 +343,17 @@
 
       ! variable 3D
-      dl_tmp3D(:,:,:)=dp_fill_sp
-
-      tg_gdept_0 =var_init('gdept_0 ',dl_tmp3D(:,:,:), dd_fill=dp_fill_sp, id_type=NF90_FLOAT)
-      tg_gdepw_0 =var_init('gdepw_0 ',dl_tmp3D(:,:,:), dd_fill=dp_fill_sp, id_type=NF90_FLOAT)
+      dl_tmp3D(:,:,:)=dp_fill
+
+      tg_gdept_0 =var_init('gdept_0 ',dl_tmp3D(:,:,:), dd_fill=dp_fill, id_type=NF90_DOUBLE)
+      tg_gdepw_0 =var_init('gdepw_0 ',dl_tmp3D(:,:,:), dd_fill=dp_fill, id_type=NF90_DOUBLE)
       !tg_gdep3w_0=var_init('gdep3w_0',dl_tmp3D(:,:,:), dd_fill=dp_fill_sp, id_type=NF90_FLOAT)
 
-      dl_tmp3D(:,:,:)=dp_fill
-      
       tg_e3t_0   =var_init('e3t_0   ',dl_tmp3D(:,:,:), dd_fill=dp_fill, id_type=NF90_DOUBLE)
       tg_e3u_0   =var_init('e3u_0   ',dl_tmp3D(:,:,:), dd_fill=dp_fill, id_type=NF90_DOUBLE)
       tg_e3v_0   =var_init('e3v_0   ',dl_tmp3D(:,:,:), dd_fill=dp_fill, id_type=NF90_DOUBLE)
       tg_e3w_0   =var_init('e3w_0   ',dl_tmp3D(:,:,:), dd_fill=dp_fill, id_type=NF90_DOUBLE)
-      !tg_e3f_0   =var_init('e3f_0   ',dl_tmp3D(:,:,:), dd_fill=dp_fill, id_type=NF90_DOUBLE)
-      !tg_e3uw_0  =var_init('e3uw_0  ',dl_tmp3D(:,:,:), dd_fill=dp_fill, id_type=NF90_DOUBLE)
-      !tg_e3vw_0  =var_init('e3vw_0  ',dl_tmp3D(:,:,:), dd_fill=dp_fill, id_type=NF90_DOUBLE)
+      tg_e3f_0   =var_init('e3f_0   ',dl_tmp3D(:,:,:), dd_fill=dp_fill, id_type=NF90_DOUBLE)
+      tg_e3uw_0  =var_init('e3uw_0  ',dl_tmp3D(:,:,:), dd_fill=dp_fill, id_type=NF90_DOUBLE)
+      tg_e3vw_0  =var_init('e3vw_0  ',dl_tmp3D(:,:,:), dd_fill=dp_fill, id_type=NF90_DOUBLE)
 
    END SUBROUTINE grid_zgr_init
@@ -409,7 +413,7 @@
       CALL var_clean(tg_e3u_0   )
       CALL var_clean(tg_e3v_0   )
-      !CALL var_clean(tg_e3f_0   )
-      !CALL var_clean(tg_e3uw_0  )
-      !CALL var_clean(tg_e3vw_0  )
+      CALL var_clean(tg_e3f_0   )
+      CALL var_clean(tg_e3uw_0  )
+      CALL var_clean(tg_e3vw_0  )
       
    END SUBROUTINE grid_zgr_clean
@@ -493,6 +497,6 @@
       REAL(dp)          :: dn_zb_b     = NF90_FILL_DOUBLE
 
-      ! namcla
-      INTEGER(i4)       :: in_cla      = 0
+!      ! namcla
+!      INTEGER(i4)       :: in_cla      = 0
 
       ! namwd
@@ -507,4 +511,5 @@
 !      LOGICAL           :: ln_zoom     = .FALSE.
       LOGICAL           :: ln_c1d      = .FALSE.
+      LOGICAL           :: ln_e3_dep   = .FALSE.
 
 !      ! namzoom
@@ -568,6 +573,6 @@
       &   dn_zb_b         !<  offset for calculating Zb
 
-      NAMELIST /namcla/ &
-      &  in_cla            !< =1 cross land advection for exchanges through some straits (ORCA2)
+!      NAMELIST /namcla/ &
+!      &  in_cla            !< =1 cross land advection for exchanges through some straits (ORCA2)
 
       NAMELIST /namwd/ &  !< wetting and drying
@@ -581,5 +586,6 @@
 !      &  in_bench,      &  !< benchmark parameter (in_mshhgr = 5 )
 !      &  ln_zoom,       &  !< use zoom
-      &  ln_c1d            !< use configuration 1D
+      &  ln_c1d,         &  !< use configuration 1D
+      &  ln_e3_dep          !< new vertical scale factors [T, F:old definition]
 
 !      NAMELIST /namzoom/&
@@ -616,5 +622,5 @@
       IF( ln_zps    ) READ( il_fileid, NML = namzps  )
       IF( ln_sco    ) READ( il_fileid, NML = namsco  )
-      READ( il_fileid, NML = namcla  )
+!      READ( il_fileid, NML = namcla  )
       READ( il_fileid, NML = namwd   )
       READ( il_fileid, NML = namgrd  )
@@ -676,5 +682,5 @@
       grid_zgr_nam%d_zb_b     = dn_zb_b
 
-      grid_zgr_nam%i_cla      = in_cla
+!      grid_zgr_nam%i_cla      = in_cla
 
       grid_zgr_nam%d_wdmin1   = dn_wdmin1  
@@ -687,4 +693,5 @@
 !      grid_zgr_nam%l_zoom     = ln_zoom
       grid_zgr_nam%l_c1d      = ln_c1d
+      grid_zgr_nam%l_e3_dep   = ln_e3_dep
 
 !      grid_zgr_nam%c_cfz      = cn_cfz   
@@ -743,4 +750,5 @@
       CALL logger_info(' s- or hybrid z-s-coordinate    ln_sco    = '//TRIM(fct_str(td_nam%l_sco)))
       CALL logger_info(' ice shelf cavities             ln_isfcav = '//TRIM(fct_str(td_nam%l_isfcav)))
+      CALL logger_info(' vertical scale factors         ln_e3_dep = '//TRIM(fct_str(td_nam%l_e3_dep)))
 
       il_count=0
@@ -759,4 +767,7 @@
       ENDIF
 
+      IF(.NOT. td_nam%l_e3_dep )THEN
+         CALL logger_info("Obsolescent definition of e3 scale factors is used")
+      ENDIF
       ! Build the vertical coordinate system
       ! ------------------------------------
@@ -773,8 +784,8 @@
       ! z-coordinate
       IF( td_nam%l_zco ) CALL grid_zgr__zco(jpk)
-      
+
       ! Partial step z-coordinate
       IF( td_nam%l_zps ) CALL grid_zgr__zps_fill( td_nam,jpi,jpj,jpk,td_bathy,td_risfdep )
-      
+
       ! s-coordinate or hybrid z-s coordinate
       IF( td_nam%l_sco ) CALL grid_zgr__sco_fill( td_nam,jpi,jpj,jpk,td_bathy )
@@ -807,9 +818,9 @@
          !&   ' 3w '//TRIM(fct_str(MINVAL( tg_gdep3w_0%d_value(:,:,:,1) )))//&
          &   '  t '//TRIM(fct_str(MINVAL( tg_e3t_0%d_value(:,:,:,1) )))//&
-         !&   '  f '//TRIM(fct_str(MINVAL( tg_e3f_0%d_value(:,:,:,1) )))//&
+         &   '  f '//TRIM(fct_str(MINVAL( tg_e3f_0%d_value(:,:,:,1) )))//&
          &   '  u '//TRIM(fct_str(MINVAL( tg_e3u_0%d_value(:,:,:,1) )))//&
          &   '  v '//TRIM(fct_str(MINVAL( tg_e3v_0%d_value(:,:,:,1) )))//&
-         !&   ' uw '//TRIM(fct_str(MINVAL( tg_e3uw_0%d_value(:,:,:,1) )))//&
-         !&   ' vw '//TRIM(fct_str(MINVAL( tg_e3vw_0%d_value(:,:,:,1) )))//&
+         &   ' uw '//TRIM(fct_str(MINVAL( tg_e3uw_0%d_value(:,:,:,1) )))//&
+         &   ' vw '//TRIM(fct_str(MINVAL( tg_e3vw_0%d_value(:,:,:,1) )))//&
          &   '  w '//TRIM(fct_str(MINVAL( tg_e3w_0%d_value(:,:,:,1) ))) )
       CALL logger_info(' MAX val depth t '//TRIM(fct_str(MAXVAL( tg_gdept_0%d_value(:,:,:,1) )))//&
@@ -817,9 +828,9 @@
          !&   ' 3w '//TRIM(fct_str(MAXVAL( tg_gdep3w_0%d_value(:,:,:,1) ))) )
       CALL logger_info(' MAX val e3    t '//TRIM(fct_str(MAXVAL( tg_e3t_0%d_value(:,:,:,1) )))//&
-         !&   ' f '//TRIM(fct_str(MAXVAL( tg_e3f_0%d_value(:,:,:,1) )))//&
+         &   ' f '//TRIM(fct_str(MAXVAL( tg_e3f_0%d_value(:,:,:,1) )))//&
          &   ' u '//TRIM(fct_str(MAXVAL( tg_e3u_0%d_value(:,:,:,1) )))//&
          &   ' v '//TRIM(fct_str(MAXVAL( tg_e3v_0%d_value(:,:,:,1) )))//&
-         !&   ' uw '//TRIM(fct_str(MAXVAL( tg_e3uw_0%d_value(:,:,:,1) )))//&
-         !&   ' vw '//TRIM(fct_str(MAXVAL( tg_e3vw_0%d_value(:,:,:,1) )))//&
+         &   ' uw '//TRIM(fct_str(MAXVAL( tg_e3uw_0%d_value(:,:,:,1) )))//&
+         &   ' vw '//TRIM(fct_str(MAXVAL( tg_e3vw_0%d_value(:,:,:,1) )))//&
          &   ' w '//TRIM(fct_str(MAXVAL( tg_e3w_0%d_value(:,:,:,1) ))) )
 
@@ -889,7 +900,7 @@
          &  td_nam%d_ppsur == dp_pp_to_be_computed           ) THEN
          !
-         za1  = (  zdzmin - zhmax / FLOAT(jpk-1)  )                                                 &
-            & / ( TANH((1-zkth)/zacr) - zacr/FLOAT(jpk-1) * (  LOG( COSH( (jpk - zkth) / zacr) )    &
-            &                                                - LOG( COSH( ( 1  - zkth) / zacr) )  )  )
+         za1  = (  zdzmin - zhmax / REAL(jpk-1,dp)  )                                               &
+            & / ( TANH((1-zkth)/zacr) - zacr/REAL(jpk-1,dp) * ( LOG( COSH( (jpk - zkth) / zacr) )   &
+            &                                               - LOG( COSH( ( 1  - zkth) / zacr) ) ) )
          za0  = zdzmin  - za1 *             TANH( (1-zkth) / zacr )
          zsur =   - za0 - za1 * zacr * LOG( COSH( (1-zkth) / zacr )  )
@@ -904,10 +915,10 @@
       CALL logger_info(' GRID ZGR Z : Reference vertical z-coordinates')
       CALL logger_info('~~~~~~~~~~~')
-      IF( zkth == 0._wp ) THEN 
+      IF( zkth == 0._dp ) THEN 
          CALL logger_info('Uniform grid with '//TRIM(fct_str(jpk-1))//' layers')
          CALL logger_info('Total depth    :'//TRIM(fct_str(zhmax)))
          CALL logger_info('Layer thickness:'//TRIM(fct_str(zhmax/(jpk-1))))         
       ELSE
-         IF( za1 == 0._wp .AND. za0 == 0._wp .AND. zsur == 0._wp ) THEN
+         IF( za1 == 0._dp .AND. za0 == 0._dp .AND. zsur == 0._dp ) THEN
             CALL logger_info('zsur, za0, za1 computed from ')
             CALL logger_info('        zdzmin = '//TRIM(fct_str(zdzmin)))
@@ -931,8 +942,8 @@
       ! init
       IF( zkth == 0._dp ) THEN            !  uniform vertical grid       
-         za1 = zhmax / FLOAT(jpk-1) 
+         za1 = zhmax / REAL(jpk-1,dp) 
          DO jk = 1, jpk
-            zw = FLOAT( jk )
-            zt = FLOAT( jk ) + 0.5_dp
+            zw = REAL( jk, dp )
+            zt = REAL( jk, dp ) + 0.5_dp
             tg_gdepw_1d%d_value(1,1,jk,1) = ( zw - 1 ) * za1
             tg_gdept_1d%d_value(1,1,jk,1) = ( zt - 1 ) * za1
@@ -943,6 +954,6 @@
          IF( .NOT. td_nam%l_dbletanh ) THEN
             DO jk = 1, jpk
-               zw = REAL( jk , wp )
-               zt = REAL( jk , wp ) + 0.5_dp
+               zw = REAL( jk , dp )
+               zt = REAL( jk , dp ) + 0.5_dp
                tg_gdepw_1d%d_value(1,1,jk,1) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) )  )
                tg_gdept_1d%d_value(1,1,jk,1) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) )  )
@@ -952,6 +963,6 @@
          ELSE
             DO jk = 1, jpk
-               zw = FLOAT( jk )
-               zt = FLOAT( jk ) + 0.5_dp
+               zw = REAL( jk, dp )
+               zt = REAL( jk, dp ) + 0.5_dp
                ! Double tanh function
                tg_gdepw_1d%d_value(1,1,jk,1) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth ) / zacr  ) )    &
@@ -968,5 +979,5 @@
       ENDIF      
 
-      IF ( td_nam%l_isfcav ) THEN
+      IF ( td_nam%l_isfcav .OR. td_nam%l_e3_dep ) THEN
          ! need to be like this to compute the pressure gradient with ISF. 
          ! If not, level beneath the ISF are not aligned (sum(e3t) /= depth)
@@ -989,5 +1000,5 @@
 !      
 !      ! ref. depth with tolerance (10% of minimum layer thickness)
-!      zrefdep = 10._wp - 0.1_wp * MINVAL( tg_e3w_1d%d_value(1,1,:,1) )
+!      zrefdep = 10._dp - 0.1_dp * MINVAL( tg_e3w_1d%d_value(1,1,:,1) )
 !      
 !      ! shallowest W level Below ~10m
@@ -1126,7 +1137,7 @@
             WHERE ( td_bathy%d_value(:,:,1,1) <= td_risfdep%d_value(:,:,1,1) + td_nam%d_isfhmin )
                tg_misfdep%d_value(:,:,1,1) = 0
-               td_risfdep%d_value(:,:,1,1) = 0._wp
+               td_risfdep%d_value(:,:,1,1) = 0._dp
                tg_mbathy%d_value (:,:,1,1) = 0
-               td_bathy%d_value  (:,:,1,1) = 0._wp
+               td_bathy%d_value  (:,:,1,1) = 0._dp
             END WHERE
          END IF
@@ -1211,8 +1222,8 @@
          tg_e3u_0%d_value   (:,:,jk,1) = tg_e3t_1d%d_value  (1,1,jk,1)
          tg_e3v_0%d_value   (:,:,jk,1) = tg_e3t_1d%d_value  (1,1,jk,1)
-         !tg_e3f_0%d_value   (:,:,jk,1) = tg_e3t_1d%d_value  (1,1,jk,1)
+         tg_e3f_0%d_value   (:,:,jk,1) = tg_e3t_1d%d_value  (1,1,jk,1)
          tg_e3w_0%d_value   (:,:,jk,1) = tg_e3w_1d%d_value  (1,1,jk,1)
-         !tg_e3uw_0%d_value  (:,:,jk,1) = tg_e3w_1d%d_value  (1,1,jk,1)
-         !tg_e3vw_0%d_value  (:,:,jk,1) = tg_e3w_1d%d_value  (1,1,jk,1)
+         tg_e3uw_0%d_value  (:,:,jk,1) = tg_e3w_1d%d_value  (1,1,jk,1)
+         tg_e3vw_0%d_value  (:,:,jk,1) = tg_e3w_1d%d_value  (1,1,jk,1)
       END DO
 
@@ -1760,27 +1771,4 @@
                      tg_gdept_0%d_value(ji,jj,ik+1,1) = tg_gdept_0%d_value(ji,jj,ik,1) + tg_e3t_0%d_value(ji,jj,ik,1)
                   ENDIF
-                  !gm Bug?  check the gdepw_1d
-                  !       ... on ik
-                  tg_gdept_0%d_value(ji,jj,ik,1) = tg_gdepw_1d%d_value ( 1, 1,ik  ,1)   &
-                                                &    + ( tg_gdepw_0%d_value  (ji,jj,ik+1,1) - tg_gdepw_1d%d_value(1,1,ik,1) )&
-                                                &    * ( (tg_gdept_1d%d_value( 1, 1,ik  ,1) - tg_gdepw_1d%d_value(1,1,ik,1)) &
-                                                &      / (tg_gdepw_1d%d_value( 1, 1,ik+1,1) - tg_gdepw_1d%d_value(1,1,ik,1)) )
-
-                  tg_e3t_0%d_value  (ji,jj,ik,1) = tg_e3t_1d%d_value  ( 1, 1,ik  ,1)   &
-                                                &    * ( tg_gdepw_0%d_value (ji,jj,ik+1,1) - tg_gdepw_1d%d_value(1,1,ik,1) ) &
-                                                &    / ( tg_gdepw_1d%d_value( 1, 1,ik+1,1) - tg_gdepw_1d%d_value(1,1,ik,1) ) 
-
-                  tg_e3w_0%d_value  (ji,jj,ik,1) = 0.5_dp   &
-                                                &     *  ( tg_gdepw_0%d_value (ji,jj,ik+1,1)   &
-                                                &        + tg_gdepw_1d%d_value( 1, 1,ik+1,1)   &
-                                                &        - tg_gdepw_1d%d_value( 1, 1,ik  ,1)*2._dp) &
-                                                &     *  ( tg_e3w_1d%d_value  ( 1, 1,ik  ,1)   &
-                                                &        / ( tg_gdepw_1d%d_value( 1, 1,ik+1,1)   &
-                                                &          - tg_gdepw_1d%d_value( 1, 1,ik  ,1) ) )
-
-                  !       ... on ik+1
-                  tg_e3w_0%d_value  (ji,jj,ik+1,1) = tg_e3t_0%d_value  (ji,jj,ik,1)
-                  tg_e3t_0%d_value  (ji,jj,ik+1,1) = tg_e3t_0%d_value  (ji,jj,ik,1)
-                  tg_gdept_0%d_value(ji,jj,ik+1,1) = tg_gdept_0%d_value(ji,jj,ik,1) + tg_e3t_0%d_value(ji,jj,ik,1)
                ENDIF
             END DO
@@ -1824,6 +1812,6 @@
          tg_e3u_0%d_value (:,:,jk,1) = tg_e3t_1d%d_value(1,1,jk,1)
          tg_e3v_0%d_value (:,:,jk,1) = tg_e3t_1d%d_value(1,1,jk,1)
-         !tg_e3uw_0%d_value(:,:,jk,1) = tg_e3w_1d%d_value(1,1,jk,1)
-         !tg_e3vw_0%d_value(:,:,jk,1) = tg_e3w_1d%d_value(1,1,jk,1)
+         tg_e3uw_0%d_value(:,:,jk,1) = tg_e3w_1d%d_value(1,1,jk,1)
+         tg_e3vw_0%d_value(:,:,jk,1) = tg_e3w_1d%d_value(1,1,jk,1)
       END DO
 
@@ -1834,20 +1822,20 @@
                tg_e3u_0%d_value (ji,jj,jk,1) = MIN( tg_e3t_0%d_value(ji,jj,jk,1), tg_e3t_0%d_value(ji+1,jj  ,jk,1) )
                tg_e3v_0%d_value (ji,jj,jk,1) = MIN( tg_e3t_0%d_value(ji,jj,jk,1), tg_e3t_0%d_value(ji  ,jj+1,jk,1) )
-               !tg_e3uw_0%d_value(ji,jj,jk,1) = MIN( tg_e3w_0%d_value(ji,jj,jk,1), tg_e3w_0%d_value(ji+1,jj  ,jk,1) )
-               !tg_e3vw_0%d_value(ji,jj,jk,1) = MIN( tg_e3w_0%d_value(ji,jj,jk,1), tg_e3w_0%d_value(ji  ,jj+1,jk,1) )
+               tg_e3uw_0%d_value(ji,jj,jk,1) = MIN( tg_e3w_0%d_value(ji,jj,jk,1), tg_e3w_0%d_value(ji+1,jj  ,jk,1) )
+               tg_e3vw_0%d_value(ji,jj,jk,1) = MIN( tg_e3w_0%d_value(ji,jj,jk,1), tg_e3w_0%d_value(ji  ,jj+1,jk,1) )
             END DO
          END DO
       END DO
 
-      !IF ( td_nam%l_isfcav ) THEN
-      !   ! (ISF) define e3uw (adapted for 2 cells in the water column)
-      !   CALL grid_zgr__isf_fill_e3uw(jpi,jpj)
-      !END IF
+      IF ( td_nam%l_isfcav ) THEN
+         ! (ISF) define e3uw (adapted for 2 cells in the water column)
+         CALL grid_zgr__isf_fill_e3uw(jpi,jpj)
+      END IF
 
       ! lateral boundary conditions
       CALL lbc_lnk( tg_e3u_0%d_value (:,:,:,1), 'U', td_nam%i_perio, 1._dp )
       CALL lbc_lnk( tg_e3v_0%d_value (:,:,:,1), 'V', td_nam%i_perio, 1._dp )
-      !CALL lbc_lnk( tg_e3uw_0%d_value(:,:,:,1), 'U', td_nam%i_perio, 1._dp )
-      !CALL lbc_lnk( tg_e3vw_0%d_value(:,:,:,1), 'V', td_nam%i_perio, 1._dp )
+      CALL lbc_lnk( tg_e3uw_0%d_value(:,:,:,1), 'U', td_nam%i_perio, 1._dp )
+      CALL lbc_lnk( tg_e3vw_0%d_value(:,:,:,1), 'V', td_nam%i_perio, 1._dp )
 
       ! set to z-scale factor if zero (i.e. along closed boundaries)
@@ -1855,6 +1843,6 @@
          WHERE( tg_e3u_0%d_value (:,:,jk,1) == 0._dp )   tg_e3u_0%d_value (:,:,jk,1) = tg_e3t_1d%d_value(1,1,jk,1)
          WHERE( tg_e3v_0%d_value (:,:,jk,1) == 0._dp )   tg_e3v_0%d_value (:,:,jk,1) = tg_e3t_1d%d_value(1,1,jk,1)
-         !WHERE( tg_e3uw_0%d_value(:,:,jk,1) == 0._dp )   tg_e3uw_0%d_value(:,:,jk,1) = tg_e3w_1d%d_value(1,1,jk,1)
-         !WHERE( tg_e3vw_0%d_value(:,:,jk,1) == 0._dp )   tg_e3vw_0%d_value(:,:,jk,1) = tg_e3w_1d%d_value(1,1,jk,1)
+         WHERE( tg_e3uw_0%d_value(:,:,jk,1) == 0._dp )   tg_e3uw_0%d_value(:,:,jk,1) = tg_e3w_1d%d_value(1,1,jk,1)
+         WHERE( tg_e3vw_0%d_value(:,:,jk,1) == 0._dp )   tg_e3vw_0%d_value(:,:,jk,1) = tg_e3w_1d%d_value(1,1,jk,1)
       END DO
 
@@ -2000,17 +1988,17 @@
 
       ! remove very shallow ice shelf (less than ~ 10m if 75L)
-      WHERE( td_risfdep%d_value(:,:,1,1) <= 10._wp .AND. &
+      WHERE( td_risfdep%d_value(:,:,1,1) <= 10._dp .AND. &
            & tg_misfdep%d_value(:,:,1,1) > 1 )
          tg_misfdep%d_value(:,:,1,1) = 0
-         td_risfdep%d_value(:,:,1,1) = 0.0_wp
+         td_risfdep%d_value(:,:,1,1) = 0.0_dp
          tg_mbathy%d_value (:,:,1,1) = 0
-         td_bathy%d_value  (:,:,1,1) = 0.0_wp
+         td_bathy%d_value  (:,:,1,1) = 0.0_dp
       END WHERE
-      WHERE( td_bathy%d_value(:,:,1,1) <= 30.0_wp .AND. &
-           & tg_gphit%d_value(:,:,1,1) < -60._wp )
+      WHERE( td_bathy%d_value(:,:,1,1) <= 30.0_dp .AND. &
+           & tg_gphit%d_value(:,:,1,1) < -60._dp )
          tg_misfdep%d_value(:,:,1,1) = 0
-         td_risfdep%d_value(:,:,1,1) = 0.0_wp
+         td_risfdep%d_value(:,:,1,1) = 0.0_dp
          tg_mbathy%d_value (:,:,1,1) = 0
-         td_bathy%d_value  (:,:,1,1) = 0.0_wp
+         td_bathy%d_value  (:,:,1,1) = 0.0_dp
       END WHERE
 
@@ -2669,6 +2657,6 @@
                   zdepwp = td_bathy%d_value(ji,jj,1,1)
                   ze3tp  = td_bathy%d_value(ji,jj,1,1) - tg_gdepw_1d%d_value(1,1,ik,1)
-                  ze3wp  = 0.5_wp * tg_e3w_1d%d_value(1,1,ik,1) &
-                     &            * ( 1._wp + ( ze3tp/tg_e3t_1d%d_value(1,1,ik,1) ) )
+                  ze3wp  = 0.5_dp * tg_e3w_1d%d_value(1,1,ik,1) &
+                     &            * ( 1._dp + ( ze3tp/tg_e3t_1d%d_value(1,1,ik,1) ) )
                   tg_e3t_0%d_value  (ji,jj,ik  ,1) = ze3tp
                   tg_e3t_0%d_value  (ji,jj,ik+1,1) = ze3tp
@@ -2895,50 +2883,50 @@
 !
 !   END SUBROUTINE grid_zgr__isf_fill_e3x
-!   !-------------------------------------------------------------------
-!   !> @brief This subroutine define e3uw 
-!   !>    (adapted for 2 cells in the water column) for ISF case 
-!   !>
-!   !> @details
-!   !>
-!   !> @author J.Paul
-!   !> @date September, 2015 - rewrite from zgr_zps
-!   !>
-!   !> @param[in] jpi
-!   !> @param[in] jpj
-!   !-------------------------------------------------------------------
-!   SUBROUTINE grid_zgr__isf_fill_e3uw(jpi,jpj) 
-!      IMPLICIT NONE
-!      ! Argument      
-!      INTEGER(i4)   , INTENT(IN   ) :: jpi
-!      INTEGER(i4)   , INTENT(IN   ) :: jpj
-!
-!      ! local variable
-!      INTEGER(i4) :: ikb, ikt
-!
-!      ! loop indices
-!      INTEGER(i4) :: ji
-!      INTEGER(i4) :: jj
-!      !----------------------------------------------------------------
-!
-!      DO jj = 2, jpj - 1 
-!         DO ji = 2, jpi - 1 
-!
-!            ikb = MAX(tg_mbathy%d_value (ji,jj,1,1), tg_mbathy%d_value (ji+1,jj,1,1))
-!            ikt = MAX(tg_misfdep%d_value(ji,jj,1,1), tg_misfdep%d_value(ji+1,jj,1,1))
-!            IF( ikb == ikt+1 )THEN
-!               tg_e3uw_0%d_value(ji,jj,ikb,1) = MIN( tg_gdept_0%d_value(ji,jj,ikb  ,1), tg_gdept_0%d_value(ji+1,jj  ,ikb  ,1) ) - &
-!               &                                MAX( tg_gdept_0%d_value(ji,jj,ikb-1,1), tg_gdept_0%d_value(ji+1,jj  ,ikb-1,1) )
-!            ENDIF
-!            
-!            ikb = MAX( tg_mbathy%d_value (ji,jj,1,1), tg_mbathy%d_value (ji,jj+1,1,1))
-!            ikt = MAX( tg_misfdep%d_value(ji,jj,1,1), tg_misfdep%d_value(ji,jj+1,1,1))
-!            IF( ikb == ikt+1 )THEN
-!               tg_e3vw_0%d_value(ji,jj,ikb,1) = MIN( tg_gdept_0%d_value(ji,jj,ikb  ,1), tg_gdept_0%d_value(ji  ,jj+1,ikb  ,1) ) - &
-!               &                                MAX( tg_gdept_0%d_value(ji,jj,ikb-1,1), tg_gdept_0%d_value(ji  ,jj+1,ikb-1,1) )
-!            ENDIF
-!         END DO
-!      END DO
-!
-!   END SUBROUTINE grid_zgr__isf_fill_e3uw
+   !-------------------------------------------------------------------
+   !> @brief This subroutine define e3uw 
+   !>    (adapted for 2 cells in the water column) for ISF case 
+   !>
+   !> @details
+   !>
+   !> @author J.Paul
+   !> @date September, 2015 - rewrite from zgr_zps
+   !>
+   !> @param[in] jpi
+   !> @param[in] jpj
+   !-------------------------------------------------------------------
+   SUBROUTINE grid_zgr__isf_fill_e3uw(jpi,jpj) 
+      IMPLICIT NONE
+      ! Argument      
+      INTEGER(i4)   , INTENT(IN   ) :: jpi
+      INTEGER(i4)   , INTENT(IN   ) :: jpj
+
+      ! local variable
+      INTEGER(i4) :: ikb, ikt
+
+      ! loop indices
+      INTEGER(i4) :: ji
+      INTEGER(i4) :: jj
+      !----------------------------------------------------------------
+
+      DO jj = 2, jpj - 1 
+         DO ji = 2, jpi - 1 
+
+            ikb = MAX(tg_mbathy%d_value (ji,jj,1,1), tg_mbathy%d_value (ji+1,jj,1,1))
+            ikt = MAX(tg_misfdep%d_value(ji,jj,1,1), tg_misfdep%d_value(ji+1,jj,1,1))
+            IF( ikb == ikt+1 )THEN
+               tg_e3uw_0%d_value(ji,jj,ikb,1) = MIN( tg_gdept_0%d_value(ji,jj,ikb  ,1), tg_gdept_0%d_value(ji+1,jj  ,ikb  ,1) ) - &
+               &                                MAX( tg_gdept_0%d_value(ji,jj,ikb-1,1), tg_gdept_0%d_value(ji+1,jj  ,ikb-1,1) )
+            ENDIF
+            
+            ikb = MAX( tg_mbathy%d_value (ji,jj,1,1), tg_mbathy%d_value (ji,jj+1,1,1))
+            ikt = MAX( tg_misfdep%d_value(ji,jj,1,1), tg_misfdep%d_value(ji,jj+1,1,1))
+            IF( ikb == ikt+1 )THEN
+               tg_e3vw_0%d_value(ji,jj,ikb,1) = MIN( tg_gdept_0%d_value(ji,jj,ikb  ,1), tg_gdept_0%d_value(ji  ,jj+1,ikb  ,1) ) - &
+               &                                MAX( tg_gdept_0%d_value(ji,jj,ikb-1,1), tg_gdept_0%d_value(ji  ,jj+1,ikb-1,1) )
+            ENDIF
+         END DO
+      END DO
+
+   END SUBROUTINE grid_zgr__isf_fill_e3uw
 !   !-------------------------------------------------------------------
 !   !> @brief This subroutine compute gdep3w_0 (vertical sum of e3w) 
@@ -3012,5 +3000,4 @@
       INTEGER(i4), INTENT(IN) :: jpi
       INTEGER(i4), INTENT(IN) :: jpj
-!      LOGICAL    , INTENT(IN) :: ld_domcfg
 
       ! local variable
@@ -3122,8 +3109,8 @@
       REAL(dp), DIMENSION(:,:), ALLOCATABLE :: zrj
       REAL(dp), DIMENSION(:,:), ALLOCATABLE :: zhbat
-      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ztmpi1
-      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ztmpi2
-      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ztmpj1
-      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ztmpj2
+      REAL(dp), DIMENSION(:,:), ALLOCATABLE :: ztmpi1
+      REAL(dp), DIMENSION(:,:), ALLOCATABLE :: ztmpi2
+      REAL(dp), DIMENSION(:,:), ALLOCATABLE :: ztmpj1
+      REAL(dp), DIMENSION(:,:), ALLOCATABLE :: ztmpj2
 
       ! loop indices
@@ -3170,5 +3157,5 @@
       ! set maximum ocean depth
       td_bathy%d_value(:,:,1,1) = MIN( td_nam%d_sbot_max, td_bathy%d_value(:,:,1,1) )
-      IF( .NOT.td_nam%l_wd ) THEN
+      IF( .NOT. td_nam%l_wd )THEN
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -3187,9 +3174,9 @@
       zenv(:,:) = td_bathy%d_value(:,:,1,1)
 
-      IF( .NOT. td_nam%l_wd ) THEN
+      IF( .NOT. td_nam%l_wd )THEN
       ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
          DO jj = 1, jpj
             DO ji = 1, jpi
-              IF( td_bathy%d_value(ji,jj,1,1) == 0._dp ) THEN
+              IF( td_bathy%d_value(ji,jj,1,1) == 0._dp )THEN
                 iip1 = MIN( ji+1, jpi )
                 ijp1 = MIN( jj+1, jpj )
@@ -3203,5 +3190,5 @@
                    &  td_bathy%d_value(iim1,ijm1,1,1) + &
                    &  td_bathy%d_value(ji  ,ijm1,1,1) + &
-                   &  td_bathy%d_value(iip1,ijp1,1,1) ) > 0._dp ) THEN
+                   &  td_bathy%d_value(iip1,ijp1,1,1) ) > 0._dp )THEN
                   zenv(ji,jj) = td_nam%d_sbot_min
                 ENDIF
@@ -3212,5 +3199,6 @@
 
       ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
-      CALL lbc_lnk( zenv(:,:), 'T', td_nam%i_perio, 1._dp) ! ?? , 'no0' )
+      ! this is only in mpp case, so here just do nothing
+      !! CALL lbc_lnk( zenv(:,:), 'T', td_nam%i_perio, 1._dp, 'no0' ) 
 
       ! smooth the bathymetry (if required)
@@ -3226,5 +3214,5 @@
       zrfact = ( 1._dp - td_nam%d_rmax ) / ( 1._dp + td_nam%d_rmax )      
 
-      ! initialise temporary evelope depth arrays
+      ! initialise temporary envelope depth arrays
       ALLOCATE(ztmpi1(jpi,jpj))
       ALLOCATE(ztmpi2(jpi,jpj))
@@ -3236,5 +3224,5 @@
       ztmpj1(:,:) = zenv(:,:)
       ztmpj2(:,:) = zenv(:,:)
-      
+
       ! initialise temporary r-value arrays
       ALLOCATE(zri(jpi,jpj))
@@ -3266,9 +3254,9 @@
                IF( (zenv(ji  ,jj) > 0._dp) .AND. &
                  & (zenv(iip1,jj) > 0._dp) )THEN
-                  zri(ji,jj) = ( zenv(iip1,jj) - zenv(ji,jj) ) / ( zenv(iip1,jj) + zenv(ji,jj) )
+                  zri(ji,jj) = ( zenv(iip1,  jj) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
                END IF
                IF( (zenv(ji,jj  ) > 0._dp) .AND. &
                  & (zenv(ji,ijp1) > 0._dp) )THEN
-                  zrj(ji,jj) = ( zenv(ji,ijp1) - zenv(ji,jj) ) / ( zenv(ji,ijp1) + zenv(ji,jj) )
+                  zrj(ji,jj) = ( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
                END IF
                IF( zri(ji,jj) >  td_nam%d_rmax ) ztmpi1(ji  ,jj  ) = zenv(iip1,jj  ) * zrfact
@@ -3293,5 +3281,6 @@
          END DO
          ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
-         CALL lbc_lnk( zenv, 'T', td_nam%i_perio, 1._dp) ! ?? , 'no0' )
+         ! this is only in mpp case, so here just do nothing
+         !!CALL lbc_lnk( zenv, 'T', td_nam%i_perio, 1._dp, 'no0' )
       ENDDO
       !     End loop     !
@@ -3345,5 +3334,5 @@
       tg_hbatv%d_value(:,:,1,1) = td_nam%d_sbot_min
       tg_hbatf%d_value(:,:,1,1) = td_nam%d_sbot_min
-      
+
       DO jj = 1, jpj-1
         DO ji = 1, jpi-1   ! NO vector opt.
@@ -3363,14 +3352,14 @@
           DO ji = 1, jpi
             IF( ABS(tg_hbatt%d_value(ji,jj,1,1)) < td_nam%d_wdmin1 )THEN
-               tg_hbatt%d_value(ji,jj,1,1) = SIGN(1._wp, tg_hbatt%d_value(ji,jj,1,1)) * td_nam%d_wdmin1
+               tg_hbatt%d_value(ji,jj,1,1) = SIGN(1._dp, tg_hbatt%d_value(ji,jj,1,1)) * td_nam%d_wdmin1
             ENDIF
             IF( ABS(tg_hbatu%d_value(ji,jj,1,1)) < td_nam%d_wdmin1 )THEN
-               tg_hbatu%d_value(ji,jj,1,1) = SIGN(1._wp, tg_hbatu%d_value(ji,jj,1,1)) * td_nam%d_wdmin1
+               tg_hbatu%d_value(ji,jj,1,1) = SIGN(1._dp, tg_hbatu%d_value(ji,jj,1,1)) * td_nam%d_wdmin1
             ENDIF
             IF( ABS(tg_hbatv%d_value(ji,jj,1,1)) < td_nam%d_wdmin1 )THEN
-               tg_hbatv%d_value(ji,jj,1,1) = SIGN(1._wp, tg_hbatv%d_value(ji,jj,1,1)) * td_nam%d_wdmin1
+               tg_hbatv%d_value(ji,jj,1,1) = SIGN(1._dp, tg_hbatv%d_value(ji,jj,1,1)) * td_nam%d_wdmin1
             ENDIF
             IF( ABS(tg_hbatf%d_value(ji,jj,1,1)) < td_nam%d_wdmin1 )THEN
-               tg_hbatf%d_value(ji,jj,1,1) = SIGN(1._wp, tg_hbatf%d_value(ji,jj,1,1)) * td_nam%d_wdmin1
+               tg_hbatf%d_value(ji,jj,1,1) = SIGN(1._dp, tg_hbatf%d_value(ji,jj,1,1)) * td_nam%d_wdmin1
             ENDIF
           END DO
@@ -3386,5 +3375,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            IF( tg_hbatu%d_value(ji,jj,1,1) == 0._dp ) THEN
+            IF( tg_hbatu%d_value(ji,jj,1,1) == 0._dp )THEN
                !No worries about the following line when l_wd == .true.
                IF( zhbat(ji,jj) == 0._dp ) tg_hbatu%d_value(ji,jj,1,1) = td_nam%d_sbot_min
@@ -3419,10 +3408,10 @@
 
 !!bug:  key_helsinki a verifer
-      IF( .NOT.td_nam%l_wd ) THEN
+      IF( .NOT.td_nam%l_wd )THEN
          dl_hift(:,:) = MIN( dl_hift(:,:), tg_hbatt%d_value(:,:,1,1) )
          dl_hifu(:,:) = MIN( dl_hifu(:,:), tg_hbatu%d_value(:,:,1,1) )
          dl_hifv(:,:) = MIN( dl_hifv(:,:), tg_hbatv%d_value(:,:,1,1) )
          dl_hiff(:,:) = MIN( dl_hiff(:,:), tg_hbatf%d_value(:,:,1,1) )
-      ENDIf
+      ENDIF
 
       CALL logger_info(' MAX val hif   t '//TRIM(fct_str(MAXVAL( dl_hift(:,:) )))//&
@@ -3464,5 +3453,5 @@
          CALL grid_zgr__sco_s_tanh( td_nam,jpi,jpj,jpk, &
          &                          dl_scosrf, &
-         &                          dl_hift, dl_hifu, dl_hifv)!, dl_hiff )
+         &                          dl_hift, dl_hifu, dl_hifv, dl_hiff )
       ENDIF
 
@@ -3477,17 +3466,17 @@
       CALL lbc_lnk( tg_e3u_0%d_value(:,:,:,1) , 'U', td_nam%i_perio, 1._dp )
       CALL lbc_lnk( tg_e3v_0%d_value(:,:,:,1) , 'V', td_nam%i_perio, 1._dp )
-      !CALL lbc_lnk( tg_e3f_0%d_value(:,:,:,1) , 'F', td_nam%i_perio, 1._dp )
+      CALL lbc_lnk( tg_e3f_0%d_value(:,:,:,1) , 'F', td_nam%i_perio, 1._dp )
       CALL lbc_lnk( tg_e3w_0%d_value(:,:,:,1) , 'W', td_nam%i_perio, 1._dp )
-      !CALL lbc_lnk( tg_e3uw_0%d_value(:,:,:,1), 'U', td_nam%i_perio, 1._dp )
-      !CALL lbc_lnk( tg_e3vw_0%d_value(:,:,:,1), 'V', td_nam%i_perio, 1._dp )
+      CALL lbc_lnk( tg_e3uw_0%d_value(:,:,:,1), 'U', td_nam%i_perio, 1._dp )
+      CALL lbc_lnk( tg_e3vw_0%d_value(:,:,:,1), 'V', td_nam%i_perio, 1._dp )
 
       IF( .NOT. td_nam%l_wd ) THEN
-         WHERE( tg_e3t_0%d_value(:,:,:,1) ==0_dp )  tg_e3t_0%d_value(:,:,:,1) = 1._dp
-         WHERE( tg_e3u_0%d_value(:,:,:,1) ==0_dp )  tg_e3u_0%d_value(:,:,:,1) = 1._dp
-         WHERE( tg_e3v_0%d_value(:,:,:,1) ==0_dp )  tg_e3v_0%d_value(:,:,:,1) = 1._dp
-         !WHERE( tg_e3f_0%d_value(:,:,:,1) ==0_dp )  tg_e3f_0%d_value(:,:,:,1) = 1._dp
-         WHERE( tg_e3w_0%d_value(:,:,:,1) ==0_dp )  tg_e3w_0%d_value(:,:,:,1) = 1._dp
-         !WHERE( tg_e3uw_0%d_value(:,:,:,1)==0_dp )  tg_e3uw_0%d_value(:,:,:,1)= 1._dp
-         !WHERE( tg_e3vw_0%d_value(:,:,:,1)==0_dp )  tg_e3vw_0%d_value(:,:,:,1)= 1._dp
+         WHERE( tg_e3t_0%d_value(:,:,:,1) == 0_dp )  tg_e3t_0%d_value(:,:,:,1) = 1._dp
+         WHERE( tg_e3u_0%d_value(:,:,:,1) == 0_dp )  tg_e3u_0%d_value(:,:,:,1) = 1._dp
+         WHERE( tg_e3v_0%d_value(:,:,:,1) == 0_dp )  tg_e3v_0%d_value(:,:,:,1) = 1._dp
+         WHERE( tg_e3f_0%d_value(:,:,:,1) == 0_dp )  tg_e3f_0%d_value(:,:,:,1) = 1._dp
+         WHERE( tg_e3w_0%d_value(:,:,:,1) == 0_dp )  tg_e3w_0%d_value(:,:,:,1) = 1._dp
+         WHERE( tg_e3uw_0%d_value(:,:,:,1)== 0_dp )  tg_e3uw_0%d_value(:,:,:,1)= 1._dp
+         WHERE( tg_e3vw_0%d_value(:,:,:,1)== 0_dp )  tg_e3vw_0%d_value(:,:,:,1)= 1._dp
       ENDIF
 
@@ -3529,9 +3518,9 @@
 
       CALL logger_info(' MIN val e3    t '//TRIM(fct_str(MINVAL( tg_e3t_0%d_value(:,:,:,1) )))//&
-         !&                           ' f '//TRIM(fct_str(MINVAL( tg_e3f_0%d_value(:,:,:,1) )))//&
+         &                           ' f '//TRIM(fct_str(MINVAL( tg_e3f_0%d_value(:,:,:,1) )))//&
          &                           ' u '//TRIM(fct_str(MINVAL( tg_e3u_0%d_value(:,:,:,1) )))//&
          &                           ' v '//TRIM(fct_str(MINVAL( tg_e3v_0%d_value(:,:,:,1) )))//&
-         !&                          ' uw '//TRIM(fct_str(MINVAL( tg_e3uw_0%d_value(:,:,:,1) )))//&
-         !&                          ' vw '//TRIM(fct_str(MINVAL( tg_e3vw_0%d_value(:,:,:,1) )))//&
+         &                          ' uw '//TRIM(fct_str(MINVAL( tg_e3uw_0%d_value(:,:,:,1) )))//&
+         &                          ' vw '//TRIM(fct_str(MINVAL( tg_e3vw_0%d_value(:,:,:,1) )))//&
          &                           ' w '//TRIM(fct_str(MINVAL( tg_e3w_0%d_value(:,:,:,1) ))) )
 
@@ -3541,9 +3530,9 @@
 
       CALL logger_info(' MAX val e3    t '//TRIM(fct_str(MAXVAL( tg_e3t_0%d_value(:,:,:,1) )))//&
-         !&                           ' f '//TRIM(fct_str(MAXVAL( tg_e3f_0%d_value(:,:,:,1) )))//&
+         &                           ' f '//TRIM(fct_str(MAXVAL( tg_e3f_0%d_value(:,:,:,1) )))//&
          &                           ' u '//TRIM(fct_str(MAXVAL( tg_e3u_0%d_value(:,:,:,1) )))//&
          &                           ' v '//TRIM(fct_str(MAXVAL( tg_e3v_0%d_value(:,:,:,1) )))//&
-         !&                          ' uw '//TRIM(fct_str(MAXVAL( tg_e3uw_0%d_value(:,:,:,1) )))//&
-         !&                          ' vw '//TRIM(fct_str(MAXVAL( tg_e3vw_0%d_value(:,:,:,1) )))//&
+         &                          ' uw '//TRIM(fct_str(MAXVAL( tg_e3uw_0%d_value(:,:,:,1) )))//&
+         &                          ' vw '//TRIM(fct_str(MAXVAL( tg_e3vw_0%d_value(:,:,:,1) )))//&
          &                           ' w '//TRIM(fct_str(MAXVAL( tg_e3w_0%d_value(:,:,:,1) ))) )
 
@@ -3554,5 +3543,5 @@
          DO jj = 1, jpj
 
-            IF( tg_hbatt%d_value(ji,jj,1,1) > 0._dp) THEN
+            IF( tg_hbatt%d_value(ji,jj,1,1) > 0._dp )THEN
                DO jk = 1, INT(tg_mbathy%d_value(ji,jj,1,1),i4)
                   ! check coordinate is monotonically increasing
@@ -3596,5 +3585,4 @@
                   ENDIF
                ENDDO
-
             ENDIF
 
@@ -3637,10 +3625,10 @@
       REAL(dp) :: zcoefw
 
-      REAL(wp) ::   ztmpu
-      REAL(wp) ::   ztmpv
-      REAL(wp) ::   ztmpf
-      REAL(wp) ::   ztmpu1
-      REAL(wp) ::   ztmpv1
-      REAL(wp) ::   ztmpf1
+      REAL(dp) ::   ztmpu
+      REAL(dp) ::   ztmpv
+      REAL(dp) ::   ztmpf
+      REAL(dp) ::   ztmpu1
+      REAL(dp) ::   ztmpv1
+      REAL(dp) ::   ztmpf1
 
       REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: z_gsigw3
@@ -3747,7 +3735,7 @@
 
                IF( td_nam%l_wd .AND. &
-                 & ( ztmpu1 < 0._wp .OR. ABS(ztmpu) < td_nam%d_wdmin1 ) )THEN
-                  z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) )
-                  z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) )
+                 & ( ztmpu1 < 0._dp .OR. ABS(ztmpu) < td_nam%d_wdmin1 ) )THEN
+                  z_esigtu3(ji,jj,jk) = 0.5_dp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) )
+                  z_esigwu3(ji,jj,jk) = 0.5_dp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) )
                ELSE
                   z_esigtu3(ji,jj,jk) = ( tg_hbatt%d_value(ji  ,jj,1,1)*z_esigt3(ji  ,jj,jk) &
@@ -3758,7 +3746,7 @@
 
                IF( td_nam%l_wd .AND. &
-                 & ( ztmpv1 < 0._wp .OR. ABS(ztmpv) < td_nam%d_wdmin1 ) )THEN
-                  z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) )
-                  z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) )
+                 & ( ztmpv1 < 0._dp .OR. ABS(ztmpv) < td_nam%d_wdmin1 ) )THEN
+                  z_esigtv3(ji,jj,jk) = 0.5_dp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) )
+                  z_esigwv3(ji,jj,jk) = 0.5_dp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) )
                ELSE
                   z_esigtv3(ji,jj,jk) = ( tg_hbatt%d_value(ji,jj  ,1,1)*z_esigt3(ji,jj  ,jk) &
@@ -3769,6 +3757,6 @@
 
                IF( td_nam%l_wd .AND. &
-                 & ( ztmpf1 < 0._wp .OR. ABS(ztmpf) < td_nam%d_wdmin1 ) )THEN
-                  z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji  ,jj  ,jk) &
+                 & ( ztmpf1 < 0._dp .OR. ABS(ztmpf) < td_nam%d_wdmin1 ) )THEN
+                  z_esigtf3(ji,jj,jk) = 0.25_dp * ( z_esigt3(ji  ,jj  ,jk) &
                                       &           + z_esigt3(ji+1,jj  ,jk) &
                                       &           + z_esigt3(ji  ,jj+1,jk) &
@@ -3787,13 +3775,13 @@
                tg_e3v_0%d_value(ji,jj,jk,1) = ( ( tg_hbatv%d_value(ji,jj,1,1) - td_nam%d_hc )*z_esigtv3(ji,jj,jk) &
                   &                            + td_nam%d_hc / REAL(jpk-1,dp) )
-               !tg_e3f_0%d_value(ji,jj,jk,1) = ( ( tg_hbatf%d_value(ji,jj,1,1) - td_nam%d_hc ) *z_esigtf3(ji,jj,jk) &
-               !   &                            + td_nam%d_hc/REAL(jpk-1,dp) )
-               !
+               tg_e3f_0%d_value(ji,jj,jk,1) = ( ( tg_hbatf%d_value(ji,jj,1,1) - td_nam%d_hc ) *z_esigtf3(ji,jj,jk) &
+                  &                            + td_nam%d_hc/REAL(jpk-1,dp) )
+               
                tg_e3w_0%d_value (ji,jj,jk,1)= ( ( tg_hbatt%d_value(ji,jj,1,1) - td_nam%d_hc )*z_esigw3 (ji,jj,jk) &
                   &                            + td_nam%d_hc / REAL(jpk-1,dp) )
-               !tg_e3uw_0%d_value(ji,jj,jk,1)= ( ( tg_hbatu%d_value(ji,jj,1,1) - td_nam%d_hc)*z_esigwu3(ji,jj,jk) &
-               !  &                             + td_nam%d_hc/REAL(jpk-1,dp) )
-               !tg_e3vw_0%d_value(ji,jj,jk,1)= ( ( tg_hbatv%d_value(ji,jj,1,1) - td_nam%d_hc)*z_esigwv3(ji,jj,jk) &
-               !  &                             + td_nam%d_hc/REAL(jpk-1,dp) )
+               tg_e3uw_0%d_value(ji,jj,jk,1)= ( ( tg_hbatu%d_value(ji,jj,1,1) - td_nam%d_hc)*z_esigwu3(ji,jj,jk) &
+                 &                             + td_nam%d_hc/REAL(jpk-1,dp) )
+               tg_e3vw_0%d_value(ji,jj,jk,1)= ( ( tg_hbatv%d_value(ji,jj,1,1) - td_nam%d_hc)*z_esigwv3(ji,jj,jk) &
+                 &                             + td_nam%d_hc/REAL(jpk-1,dp) )
             END DO
         END DO
@@ -3850,10 +3838,10 @@
       REAL(dp) ::   zzs, zzb  ! Surface and bottom cell thickness in sigma space
 
-      REAL(wp) ::   ztmpu
-      REAL(wp) ::   ztmpv
-      REAL(wp) ::   ztmpf
-      REAL(wp) ::   ztmpu1
-      REAL(wp) ::   ztmpv1
-      REAL(wp) ::   ztmpf1
+      REAL(dp) ::   ztmpu
+      REAL(dp) ::   ztmpv
+      REAL(dp) ::   ztmpf
+      REAL(dp) ::   ztmpu1
+      REAL(dp) ::   ztmpv1
+      REAL(dp) ::   ztmpf1
 
       REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: z_gsigw3
@@ -3909,5 +3897,5 @@
               
               IF( td_nam%d_efold /= 0.0_dp )THEN
-                zsmth   = tanh( (tg_hbatt%d_value(ji,jj,1,1)- td_nam%d_hc ) / td_nam%d_efold )
+                zsmth = TANH( (tg_hbatt%d_value(ji,jj,1,1)-td_nam%d_hc ) / td_nam%d_efold )
               ELSE
                 zsmth = 1.0_dp 
@@ -3932,5 +3920,5 @@
             DO jk = 1, jpk
               z_gsigw3(ji,jj,jk) =  REAL(jk-1,dp)        /REAL(jpk-1,dp)*(td_nam%d_hc/tg_hbatt%d_value(ji,jj,1,1))
-              z_gsigt3(ji,jj,jk) = (REAL(jk-1,dp)+0.5_wp)/REAL(jpk-1,dp)*(td_nam%d_hc/tg_hbatt%d_value(ji,jj,1,1))
+              z_gsigt3(ji,jj,jk) = (REAL(jk-1,dp)+0.5_dp)/REAL(jpk-1,dp)*(td_nam%d_hc/tg_hbatt%d_value(ji,jj,1,1))
             END DO
 
@@ -3951,6 +3939,6 @@
 
           DO jk = 1, jpk
-             tg_gdept_0%d_value (ji,jj,jk,1) = (dd_scosrf(ji,jj)+tg_hbatt%d_value(ji,jj,1,1))*z_gsigt3(ji,jj,jk)
-             tg_gdepw_0%d_value (ji,jj,jk,1) = (dd_scosrf(ji,jj)+tg_hbatt%d_value(ji,jj,1,1))*z_gsigw3(ji,jj,jk)
+             tg_gdept_0%d_value(ji,jj,jk,1) = (dd_scosrf(ji,jj)+tg_hbatt%d_value(ji,jj,1,1))*z_gsigt3(ji,jj,jk)
+             tg_gdepw_0%d_value(ji,jj,jk,1) = (dd_scosrf(ji,jj)+tg_hbatt%d_value(ji,jj,1,1))*z_gsigw3(ji,jj,jk)
              !tg_gdep3w_0%d_value(ji,jj,jk,1) = (dd_scosrf(ji,jj)+tg_hbatt%d_value(ji,jj,1,1))*z_gsi3w3(ji,jj,jk)
           END DO
@@ -3978,18 +3966,18 @@
 
             IF( td_nam%l_wd .AND. &
-              & ( ztmpu1 < 0._wp .OR. ABS(ztmpu) < td_nam%d_wdmin1 ) )THEN
-               z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) )
-               z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) )
+              & ( ztmpu1 < 0._dp .OR. ABS(ztmpu) < td_nam%d_wdmin1 ) )THEN
+               z_esigtu3(ji,jj,jk) = 0.5_dp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) )
+               z_esigwu3(ji,jj,jk) = 0.5_dp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) )
             ELSE
                z_esigtu3(ji,jj,jk) = ( tg_hbatt%d_value(ji  ,jj,1,1)*z_esigt3(ji  ,jj,jk) &
-                                   & + tg_hbatt%d_value(ji+1,jj,1,1)*z_esigt3(ji+1,jj,jk)) / ztmpu
+                                   & + tg_hbatt%d_value(ji+1,jj,1,1)*z_esigt3(ji+1,jj,jk) ) / ztmpu
                z_esigwu3(ji,jj,jk) = ( tg_hbatt%d_value(ji  ,jj,1,1)*z_esigw3(ji  ,jj,jk) &
-                                   & + tg_hbatt%d_value(ji+1,jj,1,1)*z_esigw3(ji+1,jj,jk)) / ztmpu
+                                   & + tg_hbatt%d_value(ji+1,jj,1,1)*z_esigw3(ji+1,jj,jk) ) / ztmpu
             ENDIF
 
             IF( td_nam%l_wd .AND. &
-              & ( ztmpv1 < 0._wp .OR. ABS(ztmpv) < td_nam%d_wdmin1 ) )THEN
-               z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) )
-               z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) )
+              & ( ztmpv1 < 0._dp .OR. ABS(ztmpv) < td_nam%d_wdmin1 ) )THEN
+               z_esigtv3(ji,jj,jk) = 0.5_dp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) )
+               z_esigwv3(ji,jj,jk) = 0.5_dp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) )
             ELSE
                z_esigtv3(ji,jj,jk) = ( tg_hbatt%d_value(ji,jj  ,1,1)*z_esigt3(ji,jj  ,jk) &
@@ -4000,6 +3988,6 @@
 
             IF( td_nam%l_wd .AND. &
-              & ( ztmpf1 < 0._wp .OR. ABS(ztmpf) < td_nam%d_wdmin1 ) )THEN
-               z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji  ,jj  ,jk) &
+              & ( ztmpf1 < 0._dp .OR. ABS(ztmpf) < td_nam%d_wdmin1 ) )THEN
+               z_esigtf3(ji,jj,jk) = 0.25_dp * ( z_esigt3(ji  ,jj  ,jk) &
                                    &           + z_esigt3(ji+1,jj  ,jk) &
                                    &           + z_esigt3(ji  ,jj+1,jk) &
@@ -4015,9 +4003,9 @@
              tg_e3u_0%d_value(ji,jj,jk,1) = (dd_scosrf(ji,jj)+tg_hbatu%d_value(ji,jj,1,1))*z_esigtu3(ji,jj,jk)
              tg_e3v_0%d_value(ji,jj,jk,1) = (dd_scosrf(ji,jj)+tg_hbatv%d_value(ji,jj,1,1))*z_esigtv3(ji,jj,jk)
-             !tg_e3f_0%d_value(ji,jj,jk,1) = (dd_scosrf(ji,jj)+tg_hbatf%d_value(ji,jj,1,1))*z_esigtf3(ji,jj,jk)
+             tg_e3f_0%d_value(ji,jj,jk,1) = (dd_scosrf(ji,jj)+tg_hbatf%d_value(ji,jj,1,1))*z_esigtf3(ji,jj,jk)
              !
-             tg_e3w_0%d_value (ji,jj,jk,1) = tg_hbatt%d_value(ji,jj,1,1)*z_esigw3 (ji,jj,jk)
-             !tg_e3uw_0%d_value(ji,jj,jk,1) = tg_hbatu%d_value(ji,jj,1,1)*z_esigwu3(ji,jj,jk)
-             !tg_e3vw_0%d_value(ji,jj,jk,1) = tg_hbatv%d_value(ji,jj,1,1)*z_esigwv3(ji,jj,jk)
+             tg_e3w_0%d_value(ji,jj,jk,1)  = tg_hbatt%d_value(ji,jj,1,1)*z_esigw3 (ji,jj,jk)
+             tg_e3uw_0%d_value(ji,jj,jk,1) = tg_hbatu%d_value(ji,jj,1,1)*z_esigwu3(ji,jj,jk)
+             tg_e3vw_0%d_value(ji,jj,jk,1) = tg_hbatv%d_value(ji,jj,1,1)*z_esigwv3(ji,jj,jk)
           END DO
 
@@ -4028,8 +4016,8 @@
       CALL lbc_lnk(tg_e3u_0 %d_value(:,:,:,1),'T', td_nam%i_perio, 1._dp)
       CALL lbc_lnk(tg_e3v_0 %d_value(:,:,:,1),'T', td_nam%i_perio, 1._dp) 
-      !CALL lbc_lnk(tg_e3f_0 %d_value(:,:,:,1),'T', td_nam%i_perio, 1._dp)
+      CALL lbc_lnk(tg_e3f_0 %d_value(:,:,:,1),'T', td_nam%i_perio, 1._dp)
       CALL lbc_lnk(tg_e3w_0 %d_value(:,:,:,1),'T', td_nam%i_perio, 1._dp)
-      !CALL lbc_lnk(tg_e3uw_0%d_value(:,:,:,1),'T', td_nam%i_perio, 1._dp) 
-      !CALL lbc_lnk(tg_e3vw_0%d_value(:,:,:,1),'T', td_nam%i_perio, 1._dp)
+      CALL lbc_lnk(tg_e3uw_0%d_value(:,:,:,1),'T', td_nam%i_perio, 1._dp) 
+      CALL lbc_lnk(tg_e3vw_0%d_value(:,:,:,1),'T', td_nam%i_perio, 1._dp)
 
       DEALLOCATE( z_gsigw3  )
@@ -4068,5 +4056,5 @@
    SUBROUTINE grid_zgr__sco_s_tanh( td_nam,jpi,jpj,jpk, &
          &                          dd_scosrf, &
-         &                          dd_hift, dd_hifu, dd_hifv)!, dd_hiff ) 
+         &                          dd_hift, dd_hifu, dd_hifv, dd_hiff ) 
       IMPLICIT NONE
       ! Argument      
@@ -4079,5 +4067,5 @@
       REAL(dp)   , DIMENSION(:,:), INTENT(IN   ) :: dd_hifu
       REAL(dp)   , DIMENSION(:,:), INTENT(IN   ) :: dd_hifv
-!      REAL(dp)   , DIMENSION(:,:), INTENT(IN   ) :: dd_hiff
+      REAL(dp)   , DIMENSION(:,:), INTENT(IN   ) :: dd_hiff
 
       ! local variable
@@ -4140,13 +4128,13 @@
               tg_e3v_0%d_value(ji,jj,jk,1) = ( (tg_hbatv%d_value(ji,jj,1 ,1) - dd_hifv(ji,jj)) &
                  &                            * tg_esigt%d_value(1 ,1 ,jk,1) + dd_hifv(ji,jj)/REAL(jpk-1,dp) )
-              !tg_e3f_0%d_value(ji,jj,jk,1) = ( (tg_hbatf%d_value(ji,jj,1 ,1) - dd_hiff(ji,jj)) &
-              !   &                            * tg_esigt%d_value(1 ,1, jk,1) + dd_hiff(ji,jj)/REAL(jpk-1,dp) )
+              tg_e3f_0%d_value(ji,jj,jk,1) = ( (tg_hbatf%d_value(ji,jj,1 ,1) - dd_hiff(ji,jj)) &
+                 &                            * tg_esigt%d_value(1 ,1, jk,1) + dd_hiff(ji,jj)/REAL(jpk-1,dp) )
                
               tg_e3w_0%d_value (ji,jj,jk,1)= ( (tg_hbatt%d_value(ji,jj,1 ,1) - dd_hift(ji,jj)) &
                  &                            * tg_esigw%d_value(1 ,1 ,jk,1) + dd_hift(ji,jj)/REAL(jpk-1,dp) )
-              !tg_e3uw_0%d_value(ji,jj,jk,1)= ( (tg_hbatu%d_value(ji,jj,1 ,1) - dd_hifu(ji,jj)) &
-              !   &                            * tg_esigw%d_value(1 ,1 ,jk,1) + dd_hifu(ji,jj)/REAL(jpk-1,dp) )
-              !tg_e3vw_0%d_value(ji,jj,jk,1)= ( (tg_hbatv%d_value(ji,jj,1 ,1) - dd_hifv(ji,jj)) &
-              !   &                            * tg_esigw%d_value(1 ,1 ,jk,1) + dd_hifv(ji,jj)/REAL(jpk-1,dp) )
+              tg_e3uw_0%d_value(ji,jj,jk,1)= ( (tg_hbatu%d_value(ji,jj,1 ,1) - dd_hifu(ji,jj)) &
+                 &                            * tg_esigw%d_value(1 ,1 ,jk,1) + dd_hifu(ji,jj)/REAL(jpk-1,dp) )
+              tg_e3vw_0%d_value(ji,jj,jk,1)= ( (tg_hbatv%d_value(ji,jj,1 ,1) - dd_hifv(ji,jj)) &
+                 &                            * tg_esigw%d_value(1 ,1 ,jk,1) + dd_hifv(ji,jj)/REAL(jpk-1,dp) )
             END DO
          END DO
@@ -4179,8 +4167,8 @@
       !!----------------------------------------------------------------------
       !
-      pf =   (   TANH( td_nam%d_theta * ( -(pk-0.5_dp) / REAL(jpk-1) + td_nam%d_thetb )  )   &
-         &     - TANH( td_nam%d_thetb * td_nam%d_theta                                )  )   &
-         & * (   COSH( td_nam%d_theta                           )                      &
-         &     + COSH( td_nam%d_theta * ( 2._dp * td_nam%d_thetb - 1._dp ) )  )              &
+      pf =   (   TANH( td_nam%d_theta * ( -(pk-0.5_dp) / REAL(jpk-1,dp) + td_nam%d_thetb )  ) &
+         &     - TANH( td_nam%d_thetb * td_nam%d_theta                                   )  ) &
+         & * (   COSH( td_nam%d_theta                                      )                  &
+         &     + COSH( td_nam%d_theta * ( 2._dp * td_nam%d_thetb - 1._dp ) )  )               &
          & / ( 2._dp * SINH( td_nam%d_theta ) )
       !
@@ -4215,8 +4203,8 @@
       !
       IF ( td_nam%d_theta == 0 ) then      ! uniform sigma
-         pf1 = - ( pk1 - 0.5_dp ) / REAL( jpk-1 )
+         pf1 = - ( pk1 - 0.5_dp ) / REAL(jpk-1,dp)
       ELSE                        ! stretched sigma
-         pf1 =   ( 1._dp - pbb ) * ( SINH( td_nam%d_theta*(-(pk1-0.5_dp)/REAL(jpk-1)) ) ) / SINH( td_nam%d_theta )              &
-            &  + pbb * (  (TANH( td_nam%d_theta*( (-(pk1-0.5_dp)/REAL(jpk-1)) + 0.5_dp) ) - TANH( 0.5_dp * td_nam%d_theta )  )  &
+         pf1 =   ( 1._dp - pbb ) * ( SINH( td_nam%d_theta*(-(pk1-0.5_dp)/REAL(jpk-1,dp)) ) ) / SINH( td_nam%d_theta )              &
+            &  + pbb * (  (TANH( td_nam%d_theta*( (-(pk1-0.5_dp)/REAL(jpk-1,dp)) + 0.5_dp) ) - TANH( 0.5_dp * td_nam%d_theta )  )  &
             &        / ( 2._dp * TANH( 0.5_dp * td_nam%d_theta ) )  )
       ENDIF
@@ -4253,6 +4241,8 @@
       TYPE(TNAMZ), INTENT(IN   ) :: td_nam
       INTEGER(i4), INTENT(IN   ) :: jpk
-      REAL(dp)   , INTENT(IN   ) :: pk1(jpk)       ! continuous "k" coordinate
-      REAL(dp)                   :: p_gamma(jpk)   ! stretched coordinate
+      REAL(dp)   , DIMENSION(:) , INTENT(IN   ) :: pk1       ! continuous "k" coordinate
+      ! function
+      REAL(dp)   , DIMENSION(jpk) :: p_gamma   ! stretched coordinate
+      ! local variable
       REAL(dp)   , INTENT(IN   ) :: pzb           ! Bottom box depth
       REAL(dp)   , INTENT(IN   ) :: pzs           ! surface box depth
@@ -4266,6 +4256,6 @@
       !!----------------------------------------------------------------------
       !
-      zn1  =  1./(jpk-1.)
-      zn2  =  1. -  zn1
+      zn1  =  1._dp / REAL(jpk-1,dp)
+      zn2  =  1._dp -  zn1
 
       za1 = (td_nam%d_alpha+2.0_dp)*zn1**(td_nam%d_alpha+1.0_dp)-(td_nam%d_alpha+1.0_dp)*zn1**(td_nam%d_alpha+2.0_dp) 
@@ -4363,5 +4353,5 @@
                   &                     + tg_gdepw_0%d_value(ji  ,jj  ,jk+1,1)&
                   &                     - tg_gdepw_0%d_value(ji  ,jj-1,jk+1,1) ) &
-                  &                 / ( tg_gdepw_0%d_value  (ji  ,jj  ,jk  ,1)&
+                  &                   / ( tg_gdepw_0%d_value(ji  ,jj  ,jk  ,1)&
                   &                     + tg_gdepw_0%d_value(ji  ,jj-1,jk  ,1)&
                   &                     - tg_gdepw_0%d_value(ji  ,jj  ,jk+1,1)&
Index: /branches/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src/iom_cdf.f90
===================================================================
--- /branches/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src/iom_cdf.f90	(revision 7232)
+++ /branches/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src/iom_cdf.f90	(revision 7233)
@@ -2103,5 +2103,5 @@
                   & 'gcost','gcosu','gcosv','gcosf', &
                   & 'gsint','gsinu','gsinv','gsinf', &
-                  & 'mbathy','misf','isf_draft',     &
+                  & 'mbathy','misf','isf_draft','isfdraft',     &
                   & 'hbatt','hbatu','hbatv','hbatf', &
                   & 'gsigt','gsigu','gsigv','gsigf', &
@@ -2212,4 +2212,5 @@
          CALL iom_cdf__check(il_status,"IOM CDF WRITE VAR DEF: ")
       ENDIF
+      CALL logger_debug("IOM CDF WRITE VAR DEF: type = "//TRIM(fct_str(tl_var%i_type)))
 
       ! remove unuseful attribute
Index: /branches/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src/mpp.f90
===================================================================
--- /branches/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src/mpp.f90	(revision 7232)
+++ /branches/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src/mpp.f90	(revision 7233)
@@ -2550,8 +2550,10 @@
    !> and the number of processors following I and J.
    !> Then the number of sea/land processors is compute with mask
-   !
+   !>
    !> @author J.Paul
    !> @date October, 2015 - Initial version
-   !
+   !> @date October, 2016
+   !> - compare index to td_lay number of proc instead of td_mpp (bug fix)
+   !>
    !> @param[in] td_mpp mpp strcuture
    !> @param[in] id_mask   sub domain mask (sea=1, land=0)
@@ -2670,5 +2672,5 @@
 
             ! east boundary
-            IF( ji == td_mpp%i_niproc )THEN
+            IF( ji == td_lay%i_niproc )THEN
                il_lei = td_lay%i_lci(ji,jj)
             ELSE
@@ -2677,5 +2679,5 @@
 
             ! north boundary
-            IF( jj == td_mpp%i_njproc )THEN
+            IF( jj == td_lay%i_njproc )THEN
                il_lej = td_lay%i_lcj(ji,jj)
             ELSE
Index: /branches/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src/phycst.f90
===================================================================
--- /branches/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src/phycst.f90	(revision 7232)
+++ /branches/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src/phycst.f90	(revision 7233)
@@ -34,10 +34,10 @@
    PUBLIC :: dp_siday   !< sideral day                        [s]
 
-   REAL(wp), PUBLIC ::   rday = 24.*60.*60.     !: day                                [s]
-   REAL(wp), PUBLIC ::   rsiyea                 !: sideral year                       [s]
-   REAL(wp), PUBLIC ::   rsiday                 !: sideral day                        [s]
+   REAL(dp), PUBLIC ::   rday = 24.*60.*60.     !: day                                [s]
+   REAL(dp), PUBLIC ::   rsiyea                 !: sideral year                       [s]
+   REAL(dp), PUBLIC ::   rsiday                 !: sideral day                        [s]
 
    REAL(dp), PARAMETER :: dp_pi = 3.14159274101257_dp
-   REAL(dp), PARAMETER :: dp_eps = EPSILON(1._dp)
+   REAL(dp), PARAMETER :: dp_eps = 0.5 * EPSILON(1._dp)
    REAL(dp), PARAMETER :: dp_rearth = 6371229._dp
    REAL(dp), PARAMETER :: dp_deg2rad = dp_pi/180.0
@@ -45,7 +45,7 @@
 
    REAL(dp), PARAMETER :: dp_day = 24.*60.*60.     
-   REAL(dp), PARAMETER :: dp_siyea = 365.25_wp * dp_day * &
-      &  2._wp * dp_pi / 6.283076_dp
-   REAL(dp), PARAMETER :: dp_siday = dp_day / ( 1._wp + dp_day / dp_siyea )
+   REAL(dp), PARAMETER :: dp_siyea = 365.25_dp * dp_day * &
+      &  2._dp * dp_pi / 6.283076_dp
+   REAL(dp), PARAMETER :: dp_siday = dp_day / ( 1._dp + dp_day / dp_siyea )
 
    REAL(dp), PARAMETER :: dp_delta=1.e-5
Index: /branches/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src/variable.f90
===================================================================
--- /branches/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src/variable.f90	(revision 7232)
+++ /branches/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src/variable.f90	(revision 7233)
@@ -1139,4 +1139,6 @@
    !> @date July, 2015
    !> - add unit factor (to change unit)
+   !> @date November, 2016
+   !> - allow to add scalar value
    !>
    !> @param[in] cd_name         variable name
@@ -1292,6 +1294,13 @@
          dl_value(1,1,1,:) = dd_value(:)
       ELSE
-         CALL logger_fatal("VAR INIT: can not add value from variable "//&
-         &  TRIM(cd_name)//". invalid dimension to be used")
+         IF( SIZE(dd_value(:)) > 1 )THEN
+            CALL logger_fatal("VAR INIT: can not add value from variable "//&
+            &  TRIM(cd_name)//". invalid dimension to be used")
+         ELSE
+            dl_value(1,1,1,1) = dd_value(1)
+            CALL logger_warn("VAR INIT: add scalar value for variable "//&
+            &  TRIM(cd_name))
+
+         ENDIF
       ENDIF
 
