diff --git a/EnKF.F90 b/EnKF.F90
index 9ff7265..b9521a4 100644
--- a/EnKF.F90
+++ b/EnKF.F90
@@ -1,3 +1,4 @@
+! File:          EnKF.F90
 !
 ! Created:       ???
 !
diff --git a/Tools/m_fixhycom_eco_metno.F90 b/Tools/m_fixhycom_eco_metno.F90
index 9373185..62c5465 100644
--- a/Tools/m_fixhycom_eco_metno.F90
+++ b/Tools/m_fixhycom_eco_metno.F90
@@ -258,47 +258,47 @@ subroutine hybgen_weno_remap(si,pi,dpi,ci,&
       return
       end subroutine hybgen_weno_remap
 
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
+! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
 
-      integer function tracr_get_incr(char2) 
-      !function returns the number of the tracer.
-      !char => integer
-      implicit none        
-      character(len=2) :: char2
+!       integer function tracr_get_incr(char2) 
+!       !function returns the number of the tracer.
+!       !char => integer
+!       implicit none        
+!       character(len=2) :: char2
       
-      tracr_get_incr=-1
-      select case (char2)
-         case ('01')
-           tracr_get_incr=1
-	 case ('02')
-           tracr_get_incr=2
-	 case ('03')
-           tracr_get_incr=3
-	 case ('04')
-           tracr_get_incr=4
-	 case ('05')
-           tracr_get_incr=5
-	 case ('06')
-           tracr_get_incr=6
-	 case ('07')
-           tracr_get_incr=7
-	 case ('08')
-           tracr_get_incr=8
-	 case ('09')
-           tracr_get_incr=9 
-	 case ('10')
-           tracr_get_incr=10
-	 case ('11')
-           tracr_get_incr=11
-	 case default
-           print *,'tracer unknown',char2
-      end select
-      return
+!       tracr_get_incr=-1
+!       select case (char2)
+!          case ('01')
+!            tracr_get_incr=1
+! 	 case ('02')
+!            tracr_get_incr=2
+! 	 case ('03')
+!            tracr_get_incr=3
+! 	 case ('04')
+!            tracr_get_incr=4
+! 	 case ('05')
+!            tracr_get_incr=5
+! 	 case ('06')
+!            tracr_get_incr=6
+! 	 case ('07')
+!            tracr_get_incr=7
+! 	 case ('08')
+!            tracr_get_incr=8
+! 	 case ('09')
+!            tracr_get_incr=9 
+! 	 case ('10')
+!            tracr_get_incr=10
+! 	 case ('11')
+!            tracr_get_incr=11
+! 	 case default
+!            print *,'tracer unknown',char2
+!       end select
+!       return
       
-      end function
+!       end function
 
 
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
+! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
       integer function compute_kisop(temp,sal,nz) 
       !function defines which layers are isopycnal 
       implicit none             
@@ -350,6 +350,77 @@ real function sig(t,s)
       
       end function
  
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
+      
+  !     real function sig_ref(k)
+  !     !function return the value of the target density for a given layer
+  !     implicit none
+  !     integer::k
+      
+  !     select case (k) 
+  !       case (1)
+	!   sig_ref=0.1
+	! case (2)
+	!   sig_ref=0.2 
+	! case (3)
+	!   sig_ref=0.3 	  
+	! case (4)
+	!   sig_ref=0.4 	   
+ 	! case (5)
+	!   sig_ref=0.5 
+  !       case (6)
+	!   sig_ref=24.05
+	! case (7)
+	!   sig_ref=24.96 
+	! case (8)
+	!   sig_ref=25.68 	  
+	! case (9)
+	!   sig_ref=26.05 	   
+ 	! case (10)
+	!   sig_ref=26.30 	  
+  !       case (11)
+	!   sig_ref=26.60
+	! case (12)
+	!   sig_ref=26.83 
+	! case (13)
+	!   sig_ref=27.03 	  
+	! case (14)
+	!   sig_ref=27.20 	   
+ 	! case (15)
+	!   sig_ref=27.33 
+  !       case (16)
+	!   sig_ref=27.46
+	! case (17)
+	!   sig_ref=27.55 
+	! case (18)
+	!   sig_ref=27.66 	  
+	! case (19)
+	!   sig_ref=27.74 	   
+ 	! case (20)
+	!   sig_ref=27.82 	  
+	! case (21)
+	!   sig_ref=27.90
+	! case (22)
+	!   sig_ref=27.97 
+	! case (23)
+	!   sig_ref=28.01 	  
+	! case (24)
+	!   sig_ref=28.04 	   
+ 	! case (25)
+	!   sig_ref=28.07 
+  !       case (26)
+	!   sig_ref=28.09
+	! case (27)
+	!   sig_ref=28.11 
+	! case (28)
+	!   sig_ref=28.13 
+       
+  !     end select    
+	       
+  !     return	 
+      
+  !     end function
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
       
       real function sig_ref(k)
@@ -359,66 +430,110 @@ real function sig_ref(k)
       
       select case (k) 
         case (1)
-	  sig_ref=0.1
-	case (2)
-	  sig_ref=0.2 
-	case (3)
-	  sig_ref=0.3 	  
-	case (4)
-	  sig_ref=0.4 	   
- 	case (5)
-	  sig_ref=0.5 
+	        sig_ref=0.1
+	      case (2)
+	        sig_ref=0.2 
+	      case (3)
+	        sig_ref=0.3 	  
+	      case (4)
+	        sig_ref=0.4 	   
+ 	      case (5)
+	        sig_ref=0.5 
         case (6)
-	  sig_ref=24.05
-	case (7)
-	  sig_ref=24.96 
-	case (8)
-	  sig_ref=25.68 	  
-	case (9)
-	  sig_ref=26.05 	   
- 	case (10)
-	  sig_ref=26.30 	  
+	        sig_ref=0.6
+	      case (7)
+	        sig_ref=0.7
+	      case (8)
+	        sig_ref=0.8 	  
+	      case (9)
+	        sig_ref=0.9 	   
+ 	      case (10)
+	        sig_ref=1.0 	  
         case (11)
-	  sig_ref=26.60
-	case (12)
-	  sig_ref=26.83 
-	case (13)
-	  sig_ref=27.03 	  
-	case (14)
-	  sig_ref=27.20 	   
- 	case (15)
-	  sig_ref=27.33 
+	        sig_ref=24.05
+	      case (12)
+	        sig_ref=25.72 
+	      case (13)
+	        sig_ref=26.40 	  
+	      case (14)
+	        sig_ref=26.90 	   
+ 	      case (15)
+	        sig_ref=27.13 
         case (16)
-	  sig_ref=27.46
-	case (17)
-	  sig_ref=27.55 
-	case (18)
-	  sig_ref=27.66 	  
-	case (19)
-	  sig_ref=27.74 	   
- 	case (20)
-	  sig_ref=27.82 	  
-	case (21)
-	  sig_ref=27.90
-	case (22)
-	  sig_ref=27.97 
-	case (23)
-	  sig_ref=28.01 	  
-	case (24)
-	  sig_ref=28.04 	   
- 	case (25)
-	  sig_ref=28.07 
+	        sig_ref=27.19
+	      case (17)
+	        sig_ref=27.25 
+	      case (18)
+	        sig_ref=27.29 	  
+	      case (19)
+	        sig_ref=27.34 	   
+ 	      case (20)
+	        sig_ref=27.42 	  
+	      case (21)
+	        sig_ref=27.50
+	      case (22)
+	        sig_ref=27.56 
+	      case (23)
+	        sig_ref=27.63 	  
+	      case (24)
+	        sig_ref=27.66 	   
+ 	      case (25)
+	        sig_ref=27.69 
         case (26)
-	  sig_ref=28.09
-	case (27)
-	  sig_ref=28.11 
-	case (28)
-	  sig_ref=28.13 
+	        sig_ref=27.72
+	      case (27)
+	        sig_ref=27.74 
+	      case (28)
+	        sig_ref=27.76
+        case (29)
+	        sig_ref=27.78
+        case (30)
+	        sig_ref=27.80
+        case (31)
+	        sig_ref=27.82
+        case (32)
+	        sig_ref=27.85
+        case (33)
+	        sig_ref=27.87
+        case (34)
+	        sig_ref=27.90
+        case (35)
+	        sig_ref=27.93
+        case (36)
+	        sig_ref=27.95
+        case (37)
+	        sig_ref=27.97
+        case (38)
+	        sig_ref=27.99
+        case (39)
+	        sig_ref=28.01
+        case (40)
+	        sig_ref=28.02
+        case (41)
+	        sig_ref=28.03
+        case (42)
+	        sig_ref=28.04
+        case (43)
+	        sig_ref=28.05
+        case (44)
+	        sig_ref=28.06
+        case (45)
+	        sig_ref=28.07
+        case (46)
+	        sig_ref=28.08
+        case (47)
+	        sig_ref=28.09
+        case (48)
+	        sig_ref=28.10
+        case (49)
+	        sig_ref=28.11
+        case (50)
+	        sig_ref=28.12
        
       end select    
 	       
       return	 
+
       
       end function
- 
 end module
diff --git a/Tools/m_put_mod_fld_nc.F90 b/Tools/m_put_mod_fld_nc.F90
index 2db8580..8e3e3d5 100644
--- a/Tools/m_put_mod_fld_nc.F90
+++ b/Tools/m_put_mod_fld_nc.F90
@@ -278,10 +278,10 @@ subroutine fix_cice(fice,hice,sss,nx,ny,ncat,restart,icerestart,yearday,Ahice,Va
                   else
                      vicen(i,j,k)=0; vsnon(i,j,k)=0; aicen(i,j,k)=0
                   endif
-               end do
+               enddo
             else
                if (Vadjust==1) then
-                  vicen(i,j,:)=0.; vsnon(i,j,:)=0
+                  vicen(i,j,:)=0.; vsnon(i,j,:)=0.
                elseif (Vadjust==2) then   ! considering the new ice generation by DA
                   do k=1,ncat
                      if (aicen(i,j,k)>0) then
diff --git a/Tools/p_fixhycom_eco.F90 b/Tools/p_fixhycom_eco.F90
index 20a1665..5d7ed1a 100644
--- a/Tools/p_fixhycom_eco.F90
+++ b/Tools/p_fixhycom_eco.F90
@@ -88,6 +88,15 @@ program fixhycom_eco
     
     real, dimension(:,:,:), allocatable::temp,sal
     integer::kisop
+
+    character(len=8) :: varname
+    character(len=100) :: headerBIOrestart
+    character(len=24) :: dummy
+    character(len=2) :: dummy2
+    integer :: vlayer, tstep, counter, cc, idummy
+    character(len=8), allocatable :: varnames(:)
+    logical :: found 
+!    real, allocatable :: sigma(:)   
 #endif
 
    icerestart=''
@@ -112,9 +121,15 @@ program fixhycom_eco
    call parse_blkdat('idm   ','integer',rdummy,idm)
    call parse_blkdat('jdm   ','integer',rdummy,jdm)
    call parse_blkdat('kdm   ','integer',rdummy,kdm)
-#if defined (ECO)
-   call parse_blkdat('ntracr','integer',rdummy,ntracr) 	  
-#endif
+! #if defined (ECO)
+!       allocate(sigma(kdm))
+!    call parse_blkdat('ntracr','integer',rdummy,ntracr) 	  
+!       do cc=1,kdm
+!          !call parse_blkdat('sigma ','real',sigma(cc),idummy)
+!          call blkinr(sigma(cc), 'sigma ','(a6," =",f10.4," m")')
+!       enddo
+!    write(*,*)sigma
+! #endif
       
    if (idm>0 .and. idm < 1e4 .and. jdm>0 .and. jdm<1e4) then
       allocate(fld  (idm,jdm))
@@ -132,17 +147,17 @@ program fixhycom_eco
       allocate(dpold(idm,jdm,kdm))
       allocate(dp   (idm,jdm,kdm))
       allocate(press(kdm+1))
-#if defined (ECO)
-      allocate(dpfor(idm,jdm,kdm))
-      allocate(cfi(kdm,ntracr,2))
-      allocate(prsf(kdm+1))
-      allocate(tracerf(idm,jdm,kdm,ntracr))
-      allocate(trcraij(kdm,ntracr))
-      allocate(lcm(kdm))
-      allocate(dpf(kdm))
-      allocate(sal(idm,jdm,kdm))
-      allocate(temp(idm,jdm,kdm))
-#endif
+!#if defined (ECO)
+!      allocate(dpfor(idm,jdm,kdm))
+!      allocate(cfi(kdm,ntracr,2))
+!      allocate(prsf(kdm+1))
+!      allocate(tracerf(idm,jdm,kdm,ntracr))
+!      allocate(trcraij(kdm,ntracr))
+!      allocate(lcm(kdm))
+!      allocate(dpf(kdm))
+!      allocate(sal(idm,jdm,kdm))
+!      allocate(temp(idm,jdm,kdm))
+!#endif
    else
       print *,'fld allocate error'
       stop '(EnKF_postprocess)'
@@ -159,15 +174,86 @@ program fixhycom_eco
       write(*,*) 'Can not find '//restart(1:fnd-1)//'.b'
       stop '(EnKF_postprocess)'
    end if
-
+    
    print *,restart(1:fnd-1)//'.b'
    newfile='fix'//restart(1:fnd-1)
 
-
    ! Get model grid
    call get_mod_grid(modlon,modlat,depths,mindx,meandx,idm,jdm)
-
+    
 #if defined (ECO)
+! This will determine the number of biogeochemistry state variables
+! and their names
+    
+    open(unit=10, file=restart(1:fnd-1)//'.b', status='old', action='read')
+    read(10, '(A100)') headerBIOrestart
+    read(10, '(A100)') headerBIOrestart
+    
+    varname = "dummyyyy"
+    do while(varname /= "dpmixl")
+       read(10, '(A8)') varname ! reads the 1st dpmixl
+    end do
+    read(10, '(A8)') varname ! reads th 2nd dpmixl
+
+    ! And hopefully at this point, the lines correspond to biogeochemistry state
+    ! variables. Check if this is the case prior to execution
+
+    ntracr = 0
+    do while(.true.)
+        read(10, '(A8,A24,I2,A2,I1)') varname,dummy,vlayer,dummy2,tstep
+        if (vlayer==0) exit
+        if (vlayer == 1 .and. tstep == 1) then
+           ntracr = ntracr + 1
+        end if
+    end do
+
+
+    ! We do the same loop, but this time we have stored "ntracr" which is the
+    ! number of biogeochemistry state variables. The following will store the
+    ! variable names
+    rewind (10)
+    allocate(varnames(ntracr))
+
+    read(10, '(A100)') headerBIOrestart
+    read(10, '(A100)') headerBIOrestart
+
+    varname = "dummyyyy"
+    do while(varname /= "dpmixl")
+       read(10, '(A8)') varname ! reads the 1st dpmixl
+    end do
+    read(10, '(A8)') varname ! reads th 2nd dpmixl
+
+    ! And hopefully at this point, the lines correspond to biogeochemistry state
+    ! variables. Check if this is the case prior to execution
+
+    counter = 0
+    do while(.true.)
+        read(10, '(A8,A24,I2,A2,I1)') varname,dummy,vlayer,dummy2,tstep
+        if (vlayer==0) exit
+        if (vlayer == 1 .and. tstep == 1) then
+           counter = counter + 1
+           varnames(counter) = varname
+        end if
+    end do
+
+
+    close(10)
+
+!#endif
+!
+!#if defined (ECO)
+      allocate(dpfor(idm,jdm,kdm))
+      allocate(cfi(kdm,ntracr,2))
+      allocate(prsf(kdm+1))
+      allocate(tracerf(idm,jdm,kdm,ntracr))
+      allocate(trcraij(kdm,ntracr))
+      allocate(lcm(kdm))
+      allocate(dpf(kdm))
+      allocate(sal(idm,jdm,kdm))
+      allocate(temp(idm,jdm,kdm))
+!#endif
+!
+!#if defined (ECO)
    !files where are stored the forecast fields!
    restfor='forecast'//cmem
    dpthin = onem*0.001
@@ -181,21 +267,21 @@ program fixhycom_eco
       dpsum=dpsum+dp(:,:,k)
 #if defined (ECO)
       !reading of forecast fields: tracers, T and S
-      call get_mod_fld_new(trim(restfor),dpfor(:,:,k),imem,'dp      ',k,1,idm,jdm,0)
+      call get_mod_fld_new(trim(restfor),dpfor(:,:,k),imem,'dp    ',k,1,idm,jdm,0)
       do ktrcr=1,ntracr
          write(ctrcr,'(i2.2)') ktrcr
-         cfld='tracer'//ctrcr
+         cfld=varnames(ktrcr) !'tracer'//ctrcr
+         print *,cfld,imem,trim(restfor)
          call get_mod_fld_new(trim(restfor),tracerf(:,:,k,ktrcr),imem,cfld,k,1,idm,jdm,0)
       enddo
       call get_mod_fld_new(trim(restfor),temp(:,:,k),imem,'temp    ',k,1,idm,jdm,0)
-      call get_mod_fld_new(trim(restfor),sal(:,:,k),imem,'saln     ',k,1,idm,jdm,0)
+      call get_mod_fld_new(trim(restfor),sal(:,:,k),imem,'saln    ',k,1,idm,jdm,0)
 #endif      
    end do
    print *,maxval(dpsum-depths*onem)
    dpold=dp
 
 
-
    ! DP correction
    do j=1,jdm
    do i=1,idm
@@ -283,6 +369,9 @@ program fixhycom_eco
    ! Loop over restart file
    rstind=1 ! Restart index
    allok=.true.
+#if defined (ECO)
+   counter = 0
+#endif
    do while ( allok)
 
       ! Get header info from restart
@@ -313,7 +402,18 @@ program fixhycom_eco
 !         end if
 
 	 call get_mod_fld_new(restart(1:fnd-1),fld(:,:),imem,cfld,vlevel,1,idm,jdm,0)
-
+#if defined (ECO) 
+         found = .false.
+
+         do cc = 1,ntracr
+            if (.not. found) then
+               if (cfld == varnames(cc)) then
+                  found = .true.
+                  counter = cc!counter + 1
+               end if
+            endif
+         enddo 
+#endif
          if (trim(cfld)=='temp') then
 
             ! need salinity as well
@@ -343,15 +443,16 @@ program fixhycom_eco
          else if (trim(cfld)=='dp') then
             fld = dp(:,:,vlevel) ! NB, one time level 
 #if defined (ECO)  
-         else if (cfld(1:6)=='tracer') then  	    
-	    !updating the file!
-	    ktrcr=tracr_get_incr(cfld(7:8))
-	    if (ktrcr==-1)then
-	      print*,'alert tracer unknow'
-	      exit
-	    endif
-            fld(:,:)= tracerf(:,:,vlevel,ktrcr)
-	    
+      !    else if (cfld(1:6)=='tracer') then  	    
+	   !  !updating the file!
+	   !  ktrcr=tracr_get_incr(cfld(7:8))
+	   !  if (ktrcr==-1)then
+	   !    print*,'alert tracer unknow'
+	   !    exit
+	   !  endif
+      !       fld(:,:)= tracerf(:,:,vlevel,ktrcr)
+         else if (found) then
+            fld(:,:)= tracerf(:,:,vlevel,counter)
 #endif	 
 	 end if ! No correction for other fields in the hycom restart file
 
diff --git a/nfw.F90 b/nfw.F90
index e9a5f1a..d1003b7 100644
--- a/nfw.F90
+++ b/nfw.F90
@@ -597,7 +597,6 @@ subroutine nfw_get_var_double2D(fname, ncid, varid, vlevel,v,nx,ny)
     end if
   end subroutine nfw_get_var_double2D
 
-
   subroutine nfw_get_var_text(fname, ncid, varid, v)
     character*(*), intent(in) :: fname
     integer, intent(in) :: ncid