Index: branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90
===================================================================
--- branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90	(revision 10724)
+++ branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90	(revision 10812)
@@ -23,4 +23,5 @@
    USE wrk_nemo        ! Memory allocation
    USE timing          ! Timing
+   USE iom
 
    IMPLICIT NONE
@@ -31,7 +32,11 @@
 
    LOGICAL , PUBLIC ::   ln_tsd_init      !: T & S data flag
+   LOGICAL , PUBLIC ::   ln_tsd_interp    !: vertical interpolation flag
    LOGICAL , PUBLIC ::   ln_tsd_tradmp    !: internal damping toward input data flag
 
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tsd   ! structure of input SST (file informations, fields read)
+   INTEGER                                 ::   jpk_init , inum_dta
+   INTEGER                                 ::   id ,linum   ! local integers
+   INTEGER                                 ::   zdim(4)
 
    !!----------------------------------------------------------------------
@@ -53,11 +58,12 @@
       LOGICAL, INTENT(in), OPTIONAL ::   ld_tradmp   ! force the initialization when tradp is used
       !
-      INTEGER ::   ios, ierr0, ierr1, ierr2, ierr3   ! local integers
+      INTEGER ::   ios, ierr0, ierr1, ierr2, ierr3, ierr4, ierr5   ! local integers
       !!
       CHARACTER(len=100)            ::   cn_dir          ! Root directory for location of ssr files
-      TYPE(FLD_N), DIMENSION( jpts) ::   slf_i           ! array of namelist informations on the fields to read
-      TYPE(FLD_N)                   ::   sn_tem, sn_sal
+      TYPE(FLD_N), DIMENSION(jpts+2)::   slf_i           ! array of namelist informations on the fields to read
+      TYPE(FLD_N)                   ::   sn_tem, sn_sal, sn_dep, sn_msk
+      
       !!
-      NAMELIST/namtsd/   ln_tsd_init, ln_tsd_tradmp, cn_dir, sn_tem, sn_sal
+      NAMELIST/namtsd/   ln_tsd_init, ln_tsd_interp, ln_tsd_tradmp, cn_dir, sn_tem, sn_sal, sn_dep, sn_msk
       !!----------------------------------------------------------------------
       !
@@ -65,5 +71,5 @@
       !
       !  Initialisation
-      ierr0 = 0  ;  ierr1 = 0  ;  ierr2 = 0  ;  ierr3 = 0
+      ierr0 = 0  ;  ierr1 = 0  ;  ierr2 = 0  ;  ierr3 = 0  ; ierr4 = 0  ;  ierr5 = 0 
       !
       REWIND( numnam_ref )              ! Namelist namtsd in reference namelist : 
@@ -84,4 +90,5 @@
          WRITE(numout,*) '   Namelist namtsd'
          WRITE(numout,*) '      Initialisation of ocean T & S with T &S input data   ln_tsd_init   = ', ln_tsd_init
+         WRITE(numout,*) '      iInterpolation of initial conditions in the vertical ln_tsd_interp = ', ln_tsd_interp
          WRITE(numout,*) '      damping of ocean T & S toward T &S input data        ln_tsd_tradmp = ', ln_tsd_tradmp
          WRITE(numout,*)
@@ -97,23 +104,50 @@
          ln_tsd_init = .FALSE.
       ENDIF
+      IF( ln_tsd_interp .AND. ln_tsd_tradmp ) THEN
+            CALL ctl_stop( 'dta_tsd_init: Tracer damping and vertical interpolation not yet configured' )   ;   RETURN
+      ENDIF
+      IF( ln_tsd_interp .AND. LEN(TRIM(sn_msk%wname)) > 0 ) THEN
+            CALL ctl_stop( 'dta_tsd_init: Using vertical interpolation and weights files not recommended' )   ;   RETURN
+      ENDIF
       !
       !                             ! allocate the arrays (if necessary)
       IF(  ln_tsd_init .OR. ln_tsd_tradmp  ) THEN
          !
-         ALLOCATE( sf_tsd(jpts), STAT=ierr0 )
+         IF( ln_tsd_interp ) THEN
+           ALLOCATE( sf_tsd(jpts+2), STAT=ierr0 ) ! to carry the addtional depth information
+         ELSE
+           ALLOCATE( sf_tsd(jpts  ), STAT=ierr0 ) 
+         ENDIF 
          IF( ierr0 > 0 ) THEN
             CALL ctl_stop( 'dta_tsd_init: unable to allocate sf_tsd structure' )   ;   RETURN
          ENDIF
          !
-                                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 )
-                                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 )
-         !
-         IF( ierr0 + ierr1 + ierr2 + ierr3 > 0 ) THEN
+         IF( ln_tsd_interp ) THEN
+            CALL iom_open ( trim(cn_dir) // trim(sn_dep%clname), inum_dta ) 
+            id = iom_varid( inum_dta, sn_dep%clvar, zdim )
+            jpk_init = zdim(3)
+            IF(lwp) WRITE(numout,*) 'Dimension of veritcal coordinate in ICs: ', jpk_init
+            CALL iom_close( inum_dta )   ! Close the input file
+            !
+                                 ALLOCATE( sf_tsd(jp_tem)%fnow(jpi,jpj,jpk_init  ) , STAT=ierr0 )
+            IF( sn_tem%ln_tint ) ALLOCATE( sf_tsd(jp_tem)%fdta(jpi,jpj,jpk_init,2) , STAT=ierr1 )
+                                 ALLOCATE( sf_tsd(jp_sal)%fnow(jpi,jpj,jpk_init  ) , STAT=ierr2 )
+            IF( sn_sal%ln_tint ) ALLOCATE( sf_tsd(jp_sal)%fdta(jpi,jpj,jpk_init,2) , STAT=ierr3 )  
+                                 ALLOCATE( sf_tsd(jp_dep)%fnow(jpi,jpj,jpk_init  ) , STAT=ierr4 )
+                                 ALLOCATE( sf_tsd(jp_msk)%fnow(jpi,jpj,jpk_init  ) , STAT=ierr5 )
+         ELSE
+                                 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 )
+                                 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 )  
+         ENDIF ! ln_tsd_interp
+
+         !
+         IF( ierr0 + ierr1 + ierr2 + ierr3 + ierr4 + ierr5 > 0 ) THEN
             CALL ctl_stop( 'dta_tsd : unable to allocate T & S data arrays' )   ;   RETURN
          ENDIF
          !                         ! fill sf_tsd with sn_tem & sn_sal and control print
          slf_i(jp_tem) = sn_tem   ;   slf_i(jp_sal) = sn_sal
+         IF( ln_tsd_interp ) slf_i(jp_dep) = sn_dep   ;   slf_i(jp_msk) = sn_msk
          CALL fld_fill( sf_tsd, slf_i, cn_dir, 'dta_tsd', 'Temperature & Salinity data', 'namtsd', no_print )
          !
@@ -143,8 +177,7 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   ptsd   ! T & S data
       !
-      INTEGER ::   ji, jj, jk, jl, jkk   ! dummy loop indicies
-      INTEGER ::   ik, il0, il1, ii0, ii1, ij0, ij1   ! local integers
+      INTEGER ::   ji, jj, jk, jl, jk_init   ! dummy loop indicies
+      INTEGER ::   ik, il0, il1, ii0, ii1, ij0, ij1        ! local integers
       REAL(wp)::   zl, zi
-      REAL(wp), POINTER, DIMENSION(:) ::  ztp, zsp   ! 1D workspace
       !!----------------------------------------------------------------------
       !
@@ -181,51 +214,47 @@
 !!gm end
       !
-      ptsd(:,:,:,jp_tem) = sf_tsd(jp_tem)%fnow(:,:,:)    ! NO mask
-      ptsd(:,:,:,jp_sal) = sf_tsd(jp_sal)%fnow(:,:,:) 
-      !
-      IF( ln_sco ) THEN                   !==   s- or mixed s-zps-coordinate   ==!
-         !
-         CALL wrk_alloc( jpk, ztp, zsp )
-         !
-         IF( kt == nit000 .AND. lwp )THEN
-            WRITE(numout,*)
-            WRITE(numout,*) 'dta_tsd: interpolates T & S data onto the s- or mixed s-z-coordinate mesh'
-         ENDIF
-         !
-         DO jj = 1, jpj                         ! vertical interpolation of T & S
-            DO ji = 1, jpi
-               DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points
+      IF( kt == nit000 .AND. lwp )THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'dta_tsd: interpolates T & S data onto current mesh'
+      ENDIF
+      !
+      IF( ln_tsd_interp ) THEN                 ! probably should use pointers in the following to make more readable
+      !
+         DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points
+            DO jj= 1, jpj
+               DO ji= 1, jpi
                   zl = gdept_0(ji,jj,jk)
-                  IF(     zl < gdept_1d(1  ) ) THEN          ! above the first level of data
-                     ztp(jk) =  ptsd(ji,jj,1    ,jp_tem)
-                     zsp(jk) =  ptsd(ji,jj,1    ,jp_sal)
-                  ELSEIF( zl > gdept_1d(jpk) ) THEN          ! below the last level of data
-                     ztp(jk) =  ptsd(ji,jj,jpkm1,jp_tem)
-                     zsp(jk) =  ptsd(ji,jj,jpkm1,jp_sal)
-                  ELSE                                      ! inbetween : vertical interpolation between jkk & jkk+1
-                     DO jkk = 1, jpkm1                                  ! when  gdept(jkk) < zl < gdept(jkk+1)
-                        IF( (zl-gdept_1d(jkk)) * (zl-gdept_1d(jkk+1)) <= 0._wp ) THEN
-                           zi = ( zl - gdept_1d(jkk) ) / (gdept_1d(jkk+1)-gdept_1d(jkk))
-                           ztp(jk) = ptsd(ji,jj,jkk,jp_tem) + ( ptsd(ji,jj,jkk+1,jp_tem) - ptsd(ji,jj,jkk,jp_tem) ) * zi 
-                           zsp(jk) = ptsd(ji,jj,jkk,jp_sal) + ( ptsd(ji,jj,jkk+1,jp_sal) - ptsd(ji,jj,jkk,jp_sal) ) * zi
+                  IF( zl < sf_tsd(jp_dep)%fnow(ji,jj,1) ) THEN                     ! above the first level of data
+                     ptsd(ji,jj,jk,jp_tem) = sf_tsd(jp_tem)%fnow(ji,jj,1) 
+                     ptsd(ji,jj,jk,jp_sal) = sf_tsd(jp_sal)%fnow(ji,jj,1)
+                  ELSEIF( zl > sf_tsd(jp_dep)%fnow(ji,jj,jpk_init) ) THEN          ! below the last level of data
+                     ptsd(ji,jj,jk,jp_tem) = sf_tsd(jp_tem)%fnow(ji,jj,jpk_init)
+                     ptsd(ji,jj,jk,jp_sal) = sf_tsd(jp_sal)%fnow(ji,jj,jpk_init)
+                  ELSE                                                             ! inbetween : vertical interpolation between jk_init & jk_init+1
+                     DO jk_init = 1, jpk_init-1                                    ! when  gdept(jk_init) < zl < gdept(jk_init+1)
+                        IF( sf_tsd(jp_msk)%fnow(ji,jj,jk_init+1) == 0 ) THEN       ! if there is no data fill down
+                           sf_tsd(jp_tem)%fnow(ji,jj,jk_init+1) = sf_tsd(jp_tem)%fnow(ji,jj,jk_init)
+                           sf_tsd(jp_sal)%fnow(ji,jj,jk_init+1) = sf_tsd(jp_sal)%fnow(ji,jj,jk_init)
+                        ENDIF
+                        IF( (zl-sf_tsd(jp_dep)%fnow(ji,jj,jk_init)) * (zl-sf_tsd(jp_dep)%fnow(ji,jj,jk_init+1)) <= 0._wp ) THEN
+                           zi = ( zl - sf_tsd(jp_dep)%fnow(ji,jj,jk_init) ) / &
+                        &       (sf_tsd(jp_dep)%fnow(ji,jj,jk_init+1)-sf_tsd(jp_dep)%fnow(ji,jj,jk_init))
+                           ptsd(ji,jj,jk,jp_tem) = sf_tsd(jp_tem)%fnow(ji,jj,jk_init) + &
+                        &                          (sf_tsd(jp_tem)%fnow(ji,jj,jk_init+1)-sf_tsd(jp_tem)%fnow(ji,jj,jk_init)) * zi
+                           ptsd(ji,jj,jk,jp_sal) = sf_tsd(jp_sal)%fnow(ji,jj,jk_init) + &
+                        &                          (sf_tsd(jp_sal)%fnow(ji,jj,jk_init+1)-sf_tsd(jp_sal)%fnow(ji,jj,jk_init)) * zi
                         ENDIF
                      END DO
                   ENDIF
-               END DO
-               DO jk = 1, jpkm1
-                  ptsd(ji,jj,jk,jp_tem) = ztp(jk) * tmask(ji,jj,jk)     ! mask required for mixed zps-s-coord
-                  ptsd(ji,jj,jk,jp_sal) = zsp(jk) * tmask(ji,jj,jk)
-               END DO
-               ptsd(ji,jj,jpk,jp_tem) = 0._wp
-               ptsd(ji,jj,jpk,jp_sal) = 0._wp
-            END DO
+               ENDDO
+            ENDDO
          END DO
-         ! 
-         CALL wrk_dealloc( jpk, ztp, zsp )
-         ! 
+         !
+         ptsd(:,:,:,jp_tem) = ptsd(:,:,:,jp_tem) *tmask(:,:,:)
+         ptsd(:,:,:,jp_sal) = ptsd(:,:,:,jp_sal) *tmask(:,:,:)
       ELSE                                !==   z- or zps- coordinate   ==!
          !                             
-         ptsd(:,:,:,jp_tem) = ptsd(:,:,:,jp_tem) * tmask(:,:,:)    ! Mask
-         ptsd(:,:,:,jp_sal) = ptsd(:,:,:,jp_sal) * tmask(:,:,:)
+         ptsd(:,:,:,jp_tem) = sf_tsd(jp_tem)%fnow(:,:,:)  * tmask(:,:,:)  ! Mask
+         ptsd(:,:,:,jp_sal) = sf_tsd(jp_sal)%fnow(:,:,:)  * tmask(:,:,:)
          !
          IF( ln_zps ) THEN                      ! zps-coordinate (partial steps) interpolation at the last ocean level
@@ -257,4 +286,6 @@
                                         DEALLOCATE( sf_tsd(jp_sal)%fnow )     ! S arrays in the structure
          IF( sf_tsd(jp_sal)%ln_tint )   DEALLOCATE( sf_tsd(jp_sal)%fdta )
+         IF( ln_tsd_interp )            DEALLOCATE( sf_tsd(jp_dep)%fnow )     ! T arrays in the structure
+         IF( ln_tsd_interp )            DEALLOCATE( sf_tsd(jp_msk)%fnow )     ! T arrays in the structure
                                         DEALLOCATE( sf_tsd              )     ! the structure itself
       ENDIF
Index: branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/par_oce.F90
===================================================================
--- branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/par_oce.F90	(revision 10724)
+++ branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/par_oce.F90	(revision 10812)
@@ -68,4 +68,6 @@
    INTEGER, PUBLIC, PARAMETER ::   jp_tem = 1    !: indice for temperature
    INTEGER, PUBLIC, PARAMETER ::   jp_sal = 2    !: indice for salinity
+   INTEGER, PUBLIC, PARAMETER ::   jp_dep = 3    !: indice for depth
+   INTEGER, PUBLIC, PARAMETER ::   jp_msk = 4    !: indice for depth
 
    !!----------------------------------------------------------------------
