!************************************************************************
! Fortran 95 OPA Nesting tools						*
!									*
!     Copyright (C) 2005 Florian Lemarié (Florian.Lemarie@imag.fr)	*
!									*
!************************************************************************
!
MODULE agrif_interpolation
  !
  USE agrif_types
  USE io_netcdf
  USE bicubic_interp
  USE bilinear_interp
  USE agrif_extrapolation
  USE agrif_readwrite 
  !
  !
  !************************************************************************
  ! 									*
  ! MODULE  AGRIF_INTERPOLATION						*
  !									*
  ! module containing subroutine used for : 				*
  !   - Forcing data interpolation					*
  !   - Parent to Child coordinates interpolation			*
  !									*
  !************************************************************************
  !          
  ! 
CONTAINS
  !  
  !
  !****************************************************************
  !   subroutine agrif_interp					*
  !								*
  ! subroutine to interpolate coordinates		 	*
  !								*
  ! - input :							*
  !     tabin : coarse grid coordinate variable 		*
  !   typevar : position of interpolated variable on cells	*
  !								*
  ! - ouput :							*
  !    tabout : coordinate variable interpolated on fine grid	*
  !								*
  !****************************************************************
  !     
  SUBROUTINE agrif_interp(tabin,tabout,typevar)
    !      
    REAL*8,DIMENSION(:,:) :: tabin,tabout
    CHARACTER(*) :: typevar
    REAL*8 :: cff1
    INTEGER :: ii,jj
    REAL*8,DIMENSION(:,:),ALLOCATABLE :: tabouttemp
    !
    ! 
    !
    IF( ln_agrif_domain ) THEN
       CALL agrif_base_interp2(tabin,tabout,imin,jmin,typevar)
    ELSE
       CALL agrif_base_interp3(tabin,tabout,imin,jmin,typevar)
    ENDIF
    
    !
  END SUBROUTINE agrif_interp
  !
  !      
  SUBROUTINE agrif_base_interp2(tabin,tabout,i_min,j_min,typevar)
    !
    IMPLICIT NONE
    REAL*8,DIMENSION(:,:) :: tabin,tabout
    REAL*8 :: cff1
    INTEGER :: i_min,j_min
    INTEGER :: ii,jj,i,j
    CHARACTER(*) :: typevar   
    REAL*8,DIMENSION(:),ALLOCATABLE :: xc,yc,xf,yf
    REAL*8,DIMENSION(:,:),ALLOCATABLE :: tabtemp
    REAL*8 :: dxc,dyc,dxf,dyf
    REAL*8 :: decalxc,decalyc,decalxf,decalyf

    INTEGER :: ptx,pty
    REAL*8 :: val(4),xval(4)

    INTEGER :: nxc,nyc,nxf,nyf
    REAL :: xmin,ymin,x,y
    INTEGER :: ic,jc,itemp,jtemp

    nxc = SIZE(tabin,1)
    nyc = SIZE(tabin,2)

    nxf = SIZE(tabout,1)
    nyf = SIZE(tabout,2)      
    !

    ALLOCATE(xc(nxc))
    ALLOCATE(yc(nyc))
    ALLOCATE(xf(nxf))
    ALLOCATE(yf(nyf))

    ALLOCATE(tabtemp(nxc,nyf))

    dxc = 1.
    dyc = 1.
    dxf = 1./REAL(irafx)
    dyf = 1./REAL(irafy)

    IF (typevar .EQ. 'F') THEN
       ptx = 1 + nbghostcellsfine 
       pty = 1 + nbghostcellsfine 
       decalxc = 0.
       decalyc = 0.
       decalxf = 0.
       decalyf = 0.
    ELSEIF (typevar .EQ. 'T') THEN
       ptx = 2 + nbghostcellsfine 
       pty = 2 + nbghostcellsfine 
       decalxc = dxc/2.
       decalyc = dyc/2.
       decalxf = dxf/2.
       decalyf = dyf/2.
    ELSEIF (typevar .EQ. 'U') THEN
       ptx = 1 + nbghostcellsfine 
       pty = 2 + nbghostcellsfine 
       decalxc = 0.
       decalyc = dyc/2.
       decalxf = 0.
       decalyf = dyf/2.
    ELSEIF (typevar .EQ. 'V') THEN
       ptx = 2 + nbghostcellsfine 
       pty = 1 + nbghostcellsfine 
       decalxc = dxc/2.
       decalyc = 0.
       decalxf = dxf/2.
       decalyf = 0.
    ENDIF

    DO i=1,nxc
       xc(i) = 0. + (i-ptx) * dxc + decalxc
    ENDDO

    DO j=1,nyc
       yc(j) = 0. + (j-pty) * dyc + decalyc
    ENDDO


    xmin = (i_min - 1) * dxc
    ymin = (j_min - 1) * dyc

    DO i=1,nxf
       xf(i) = xmin + (i-ptx) * dxf + decalxf
    ENDDO

    DO j=1,nyf
       yf(j) = ymin + (j-pty) * dyf + decalyf
    ENDDO


    DO j = 1,nyf
       DO i = 1,nxc
	  x = xc(i)
	  y = yf(j)

	  ic = ptx + NINT((x-0.-decalxc)/dxc)
	  jc = pty + agrif_int((y-0.-decalyc)/dyc)

	  jc = jc - 1

	  CALL polint(yc(jc:jc+3),tabin(ic,jc:jc+3),4,y,tabtemp(i,j))
       ENDDO
    ENDDO

    DO j = 1,nyf
       DO i = 1,nxf
	  x = xf(i)
	  y = yf(j)

	  itemp = ptx + agrif_int((x-0.-decalxc)/dxc)
	  jtemp = pty + NINT((y-ymin-decalyf)/dyf)

	  itemp = itemp - 1

	  val = tabtemp(itemp:itemp+3,jtemp)
	  xval = xc(itemp:itemp+3)
	  CALL polint(xval,val,4,x,tabout(i,j))

       ENDDO
    ENDDO

    DEALLOCATE(xc,yc,xf,yf,tabtemp)


  END SUBROUTINE agrif_base_interp2
  !
  !      
  SUBROUTINE agrif_base_interp3(tabin,tabout,i_min,j_min,typevar)
    !
    IMPLICIT NONE
    REAL*8,DIMENSION(:,:) :: tabin,tabout
    INTEGER :: i_min,j_min
    CHARACTER(*) :: typevar   

    INTEGER :: nxf,nyf,nxc,nyc,zx,zy
    INTEGER :: ji,jj,jif,jjf,jic,jjc,jic1,jjc1,jdecx,jdecy
    REAL*8 :: Ax, Bx, Ay, By

    nxf = SIZE(tabout,1)
    nyf = SIZE(tabout,2)

    nxc = SIZE(tabin,1)
    nyc = SIZE(tabin,2)

    SELECT CASE(typevar)
    CASE('T')
       IF(MOD(irafx,2)==1) THEN ! odd
          zx = 1 ; zy = 1 ; jdecx = FLOOR(REAL(irafx)/2.) ; jdecy = FLOOR(REAL(irafy)/2.)
       ELSE                     ! even
          zx = 2 ; zy = 2 ; jdecx = FLOOR(REAL(irafx)/2.) ; jdecy = FLOOR(REAL(irafy)/2.)
       ENDIF
    CASE('U')
       IF(MOD(irafx,2)==1) THEN ! odd
          zx = 1 ; zy = 1 ; jdecx = irafx - 1 ; jdecy = FLOOR(REAL(irafy)/2.)
       ELSE                     ! even
          zx = 1 ; zy = 2 ; jdecx = irafx - 1 ; jdecy = FLOOR(REAL(irafy)/2.)
       ENDIF
    CASE('V')
       IF(MOD(irafx,2)==1) THEN ! odd
          zx = 1 ; zy = 1 ; jdecx = FLOOR(REAL(irafx)/2.) ; jdecy = irafy - 1
       ELSE                     ! even
          zx = 2 ; zy = 1 ; jdecx = FLOOR(REAL(irafx)/2.) ; jdecy = irafy - 1
       ENDIF
    CASE('F')
       IF(MOD(irafx,2)==1) THEN ! odd
          zx = 1 ; zy = 1 ; jdecx = irafx - 1 ; jdecy = irafy - 1
       ELSE                     ! even
          zx = 1 ; zy = 1 ; jdecx = irafx - 1 ; jdecy = irafy - 1
       ENDIF
    END SELECT

    
    DO jj = 1, nyf

       jjf = jj - jdecy
       jjc = j_min + FLOOR(REAL(jjf-1) / REAL(irafy)) 

       DO ji = 1, nxf
          
          jif = ji - jdecx 
          jic = i_min + FLOOR(REAL(jif-1) / REAL(irafx)) 

          Bx = REAL(MODULO( zx*jif-1, zx*irafx )) / REAL(zx*irafx)
          By = REAL(MODULO( zy*jjf-1, zy*irafy )) / REAL(zy*irafy)
          Ax = 1. - Bx
          Ay = 1. - By

          jic1 = MIN( nxc, jic+1 ) ! avoid out of bounds for tabin below
          jjc1 = MIN( nyc, jjc+1 ) !             --

          jic = MAX( 1, jic ) ! avoid out of bounds
          jjc = MAX( 1, jjc ) ! avoid out of bounds

          tabout(ji,jj) = ( Bx * tabin(jic1,jjc ) + Ax * tabin(jic,jjc ) ) * Ay + &
             &            ( Bx * tabin(jic1,jjc1) + Ax * tabin(jic,jjc1) ) * By
       END DO
    END DO
    
    !
  END SUBROUTINE agrif_base_interp3

  !       
  SUBROUTINE polint(xin,valin,n,x,val)
    IMPLICIT NONE
    INTEGER n
    REAL*8 xin(n), valin(n)
    REAL*8 x,val

    INTEGER ns,i,m
    REAL *8 dif,dift
    REAL*8 c(n),d(n),ho,hp,w,den,dy

    ns = 1
    dif = ABS(x-xin(1))

    DO i=1,n
       dift = ABS(x-xin(i))
       IF (dift < dif) THEN
	  ns = i
	  dif = dift
       ENDIF
       c(i) = valin(i)
       d(i) = valin(i)
    ENDDO

    val = valin(ns)
    ns = ns - 1

    DO m=1,n-1
       DO i=1,n-m
	  ho = xin(i)-x
	  hp = xin(i+m)-x
	  w=c(i+1)-d(i)
	  den = w/(ho-hp)
          d(i) = hp * den
	  c(i) = ho * den
       ENDDO
       IF (2*ns < (n-m)) THEN
          dy = c(ns+1)
       ELSE
          dy = d(ns)
          ns = ns - 1
       ENDIF
       val = val + dy
    ENDDO
  END SUBROUTINE polint
  !
  !      
  SUBROUTINE polcoe(xin,valin,n,cof)
    IMPLICIT NONE
    INTEGER n
    REAL*8 xin(n),valin(n),cof(n)


    INTEGER i,j,k
    REAL*8 b,ff,phi,s(n)

    s = 0.
    cof = 0.

    s(n)=-xin(1)
    DO i=2,n
       DO j=n+1-i,n-1
          s(j) =s(j) -xin(i)*s(j+1)
       ENDDO
       s(n)=s(n)-xin(i)
    ENDDO

    DO j=1,n
       phi=n
       DO k=n-1,1,-1
          phi = k*s(k+1)+xin(j)*phi
       ENDDO
       ff=valin(j)/phi
       b=1.
       DO k=n,1,-1
	  cof(k)=cof(k)+b*ff
	  b=s(k)+xin(j)*b
       ENDDO
    ENDDO

    RETURN
  END SUBROUTINE polcoe
  !

  !****************************************************************
  !   subroutine Correctforconservation				*
  !								*
  ! Conservation on domain boundaries ( over 2 coarse grid cells)	*
  !								*
  !****************************************************************          
  !
  !
  SUBROUTINE Correctforconservation(tabcoarse,tabfine,e1parent,e2parent,e1,e2,nxfin,nyfin,posvar,i_min,j_min)
    IMPLICIT NONE
    INTEGER :: nxfin,nyfin
    REAL*8,DIMENSION(:,:) :: tabcoarse,tabfine,e1parent,e2parent,e1,e2
    CHARACTER*1  :: posvar
    INTEGER :: i_min,j_min,diff
    INTEGER ji,jj,ipt,jpt,i,j
    INTEGER ind1,ind2,ind3,ind4,ind5,ind6
    REAL*8 cff1,cff2,cff3      
    !
    diff = 0      
    IF ( MOD(irafx,2) .EQ. 0 ) diff = 1
    !       
    ind1 = 2 + CEILING(irafx/2.0) + diff
    ind2 = nxfin-(2-1)-CEILING(irafx/2.0)
    ind3 = 2 + CEILING(irafy/2.0) + diff
    ind4 = nyfin-(2-1)-CEILING(irafy/2.0)
    ind5 = nxfin - 1 - irafx - CEILING(irafx/2.0) 
    ind6 = nyfin - 1 - irafy - CEILING(irafy/2.0) 
    !       
    IF (posvar.EQ.'T') THEN
       !       
       PRINT *,'correction'
       !
       DO ji=ind1,ind1+irafx,irafx
          !       
          ipt=i_min+(3-1)+(ji-ind1)/irafx
          !            
          DO jj=ind3,ind4,irafy
             !            
             jpt=j_min+(3-1)+(jj-ind3)/irafy

             cff1=SUM(e1(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)* &
                  e2(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff))

             cff2=SUM(e1(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)*   &
                  e2(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)*   &
                  tabfine(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff))

             cff3=e1parent(ipt,jpt)*e2parent(ipt,jpt)*tabcoarse(ipt,jpt)                        

             tabfine(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)=  &
                  tabfine(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)+(cff3-cff2)/cff1

          END DO

       ENDDO
       !****      
       DO ji=ind5,ind5+irafx,irafx
          !       
          ipt=i_min+(3-1)+(ji-ind1)/irafx

          DO jj=ind3,ind4,irafy
             jpt=j_min+(3-1)+(jj-ind3)/irafy

             cff1=SUM(e1(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)* &
                  e2(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff))

             cff2=SUM(e1(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)*   &
                  e2(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)*   &
                  tabfine(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff))

             cff3=e1parent(ipt,jpt)*e2parent(ipt,jpt)*tabcoarse(ipt,jpt)                        

             tabfine(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)=  &
                  tabfine(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)+(cff3-cff2)/cff1

          END DO
       ENDDO

       DO jj=ind3,ind3+irafy,irafy
          !       
          jpt=j_min+(3-1)+(jj-ind3)/irafy
          !
          DO ji=ind1,ind2,irafx
             !          
             ipt=i_min+(3-1)+(ji-ind1)/irafx

             cff1=SUM(e1(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)* &
                  e2(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff))

             cff2=SUM(e1(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)*   &
                  e2(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)*   &
                  tabfine(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff))

             cff3=e1parent(ipt,jpt)*e2parent(ipt,jpt)*tabcoarse(ipt,jpt)                        

             tabfine(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)=  &
                  tabfine(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)+(cff3-cff2)/cff1

          END DO
       ENDDO
       !****      
       DO jj=ind6,ind6+irafy,irafy
          !
          jpt=j_min+(3-1)+(jj-ind3)/irafy
          !
          DO ji=ind1,ind2,irafx          
             ipt=i_min+(3-1)+(ji-ind1)/irafx	

             cff1=SUM(e1(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)* &
                  e2(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff))

             cff2=SUM(e1(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)*   &
                  e2(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)*   &
                  tabfine(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff))

             cff3=e1parent(ipt,jpt)*e2parent(ipt,jpt)*tabcoarse(ipt,jpt)                        

             tabfine(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)=  &
                  tabfine(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)+(cff3-cff2)/cff1

          END DO
       ENDDO
       !                                     
       !       
    ENDIF
    RETURN
  END SUBROUTINE Correctforconservation

  !
  !****************************************************************
  !   subroutine Interp_Extrap_var				*
  !								*
  ! Interpolation and Extrapolation for temperature and salinity	*
  !								*
  ! -input : 							*
  !       filename : file containing variable to interpolate	*
  !								*
  !****************************************************************          
  !
  SUBROUTINE Interp_Extrap_var(filename, cd_type)      
    !
    IMPLICIT NONE
    !
    ! variables declaration
    !      
    CHARACTER(len=80),INTENT(in) :: filename
    CHARACTER(len=1 ), OPTIONAL, INTENT(in) :: cd_type

    CHARACTER*100 :: interp_type
    CHARACTER*100 :: Child_file,Childcoordinates
    CHARACTER*100 :: varname,childbathy
    CHARACTER*1  :: posvar
    CHARACTER*20,DIMENSION(:),POINTER :: Ncdf_varname
    CHARACTER(len=20),DIMENSION(4) :: dimnames
    !
    REAL*8, POINTER, DIMENSION(:,:) :: lonParent,latParent => NULL()
    REAL*8, POINTER, DIMENSION(:,:) :: lonChild,latChild,latlon_temp => NULL()
    REAL*8, POINTER, DIMENSION(:,:,:,:) :: tabinterp4d,tabvar1,tabvar2,tabvar3 => NULL()
    REAL*8, POINTER, DIMENSION(:,:,:) :: tabinterp3d,tabvar3d,tabvar01,tabvar02,tabvar03 => NULL()
    REAL*8, POINTER, DIMENSION(:) :: timedepth_temp,depth => NULL()
    REAL*8,DIMENSION(:,:),POINTER :: matrix => NULL()
    INTEGER,DIMENSION(:),POINTER :: src_add,dst_add  => NULL()
    INTEGER, POINTER, DIMENSION(:) :: tabtemp1DInt => NULL()
    REAL*8, POINTER, DIMENSION(:) :: nav_lev => NULL()
    !      
    LOGICAL, DIMENSION(:,:,:), POINTER :: detected_pts => NULL()  
    REAL*8,DIMENSION(:,:,:),POINTER :: mask => NULL()
    LOGICAL,DIMENSION(:,:),POINTER :: masksrc => NULL()
    LOGICAL :: Interpolation,conservation,Pacifique,Extrapolation,land_level
    !      
    INTEGER :: deptht,time,i,status,ncid,t,ii,j,nb,numlon,numlat,unlimdimid

    !        
    TYPE(Coordinates) :: G0,G1      
    !
    status = nf90_open(TRIM(filename),NF90_NOWRITE,ncid)
    IF (status/=nf90_noerr) THEN
       WRITE(*,*)"unable to open netcdf file : ",TRIM(filename)
       STOP
    ENDIF
    !
    !*****************
    !If coarse grid is masked possibility to activate an extrapolation process
    !
    Extrapolation = .FALSE.
    PRINT*,'DEBUT INTERP_EXTRAP_VAR'
    !
    !*****************
    !
    ! check grid position
    !
    IF( PRESENT(cd_type) ) THEN
       posvar = cd_type
    ELSE
       posvar = 'T'
    ENDIF     

    !
    ! read dimensions in netcdf file
    !
    CALL Read_Ncdf_dim('x',filename,numlon)
    CALL Read_Ncdf_dim('y',filename,numlat)
    IF ( Dims_Existence( 'deptht' , filename ) ) THEN
       CALL Read_Ncdf_dim('deptht',filename,deptht)
    ELSE IF ( Dims_Existence( 'depthu' , filename ) ) THEN
       CALL Read_Ncdf_dim('depthu',filename,deptht)
    ELSE IF ( Dims_Existence( 'depthv' , filename ) ) THEN
       CALL Read_Ncdf_dim('depthv',filename,deptht)
    ELSE IF ( Dims_Existence( 'depthw' , filename ) ) THEN
       CALL Read_Ncdf_dim('depthw',filename,deptht)
    ELSE IF ( Dims_Existence( 'z' , filename ) ) THEN
       CALL Read_Ncdf_dim('z',filename,deptht)
    ELSE
       deptht = N
    ENDIF
    IF    ( Dims_Existence( 'time_counter' , filename ) ) THEN
       CALL Read_Ncdf_dim('time_counter',filename,time)
       ! check that time_counter is the unlimited dim
       status = nf90_inquire(ncid, unlimiteddimid = unlimdimid)
       IF ( unlimdimid == -1 ) THEN
          WRITE(*,*)"time_counter should be the unlimited dimension in this file : ",filename
          WRITE(*,*)" use nco command: ncks -O --mk_rec_dmn time_counter ",filename," ",filename
          STOP
       ENDIF
    ELSEIF( Dims_Existence( 'time' , filename ) ) THEN
       CALL Read_Ncdf_dim('time',filename,time)
       ! check that time is the unlimited dim
       status = nf90_inquire(ncid, unlimiteddimid = unlimdimid)
       IF ( unlimdimid == -1 ) THEN
          WRITE(*,*)"time should be the unlimited dimension in this file : ",filename
          WRITE(*,*)" use nco command: ncks -O --mk_rec_dmn time ",filename," ",filename
          STOP
       ENDIF
    ELSE
       time = 0
    ENDIF
    
    !
    ! retrieve netcdf variable name
    !
    CALL Read_Ncdf_VarName(filename,Ncdf_varname)
    !
    ! define name of child grid file
    !
    CALL set_child_name(filename,Child_file) 
    CALL set_child_name(parent_coordinate_file,Childcoordinates)
    CALL set_child_name(parent_bathy_level,childbathy)
    WRITE(*,*) 'Child grid file name = ',TRIM(Child_file)
    !
    ! create this file
    !
    status = nf90_create(Child_file,NF90_WRITE,ncid)
    status = nf90_close(ncid)
    !
    ! read coordinates of both domains
    !           
    status = Read_Coordinates(parent_coordinate_file,G0)
    status = Read_Coordinates(Childcoordinates,G1,Pacifique)
    !
    ! check consistency of informations read in namelist 
    !
    IF( imax > SIZE(G0%glamt,1) .OR. jmax > SIZE(G0%glamt,2) .OR. &
         imax <= imin .OR. jmax <= jmin ) THEN                    
       WRITE(*,*) 'ERROR ***** bad child grid definition ...'
       WRITE(*,*) 'please check imin,jmin,imax,jmax,jpizoom,jpjzoom values'       
       WRITE(*,*) ' '
       STOP
    ENDIF
    !
    IF( SIZE(G1%nav_lon,1) .NE. nxfin .OR. SIZE(G1%nav_lon,2) .NE. nyfin ) THEN
       WRITE(*,*) 'ERROR ***** bad child coordinates or bathy file ...'
       WRITE(*,*) ' '
       WRITE(*,*) 'please check that child coordinates file and child bathymetry file'
       WRITE(*,*) 'has been created with the current namelist '
       WRITE(*,*) ' '
       STOP
    ENDIF
    !      
    !
    ! Initialization of T-mask thanks to bathymetry 
    !
    IF( Extrapolation ) THEN

       WRITE(*,*) 'mask initialisation on coarse and fine grids'
       ! 
       ALLOCATE(mask(numlon,numlat,N))
       CALL Init_mask(childbathy,G1,1,1)
       CALL Init_mask(parent_bathy_level,G0,1,1)
       !
    ENDIF
    !      
    ! select coordinates to use according to variable position
    !
    ALLOCATE(lonParent(numlon,numlat),latParent(numlon,numlat))
    ALLOCATE(lonChild(nxfin,nyfin),latChild(nxfin,nyfin))

    SELECT CASE(posvar)
    CASE('T')
       IF ( Vars_Existence( 'glamt' , parent_coordinate_file ) ) THEN
          lonParent = G0%glamt
          latParent = G0%gphit
       ELSEIF( Vars_Existence( 'nav_lon' , parent_coordinate_file ) ) THEN
          lonParent = G0%nav_lon
          latParent = G0%nav_lat
       ENDIF
       IF ( Vars_Existence( 'glamt' , Childcoordinates ) ) THEN
          lonChild  = G1%glamt
          latChild  = G1%gphit
       ELSEIF( Vars_Existence( 'nav_lon' , Childcoordinates ) ) THEN
          lonChild  = G1%nav_lon
          latChild  = G1%nav_lat
       ENDIF
       !!IF ( Vars_Existence( 'tmask' , Childcoordinates ) ) THEN
          mask = G1%tmask
       !!ELSEIF( Vars_Existence( 'top_level' , Childcoordinates ) ) THEN
       !!   mask = MAX( 1., G1%top_level )
       !!ELSE
       !!   mask = 1.
       !!ENDIF
    CASE('U')
       IF ( Vars_Existence( 'glamu' , parent_coordinate_file ) ) THEN
          lonParent = G0%glamu
          latParent = G0%gphiu
       ELSE
          WRITE(*,*) 'ERROR ***** missing glamu/gphiu coordinates in ', TRIM(parent_coordinate_file)
          STOP
       ENDIF
       IF ( Vars_Existence( 'glamu' , Childcoordinates ) ) THEN
          lonChild  = G1%glamu
          latChild  = G1%gphiu
       ELSE
          WRITE(*,*) 'ERROR ***** missing glamu/gphiu coordinates in ', TRIM(Childcoordinates)
          STOP
       ENDIF
       !!IF ( Vars_Existence( 'umask' , Childcoordinates ) ) THEN
          mask = G1%umask
       !!ELSE
       !!   mask = 1.
       !!ENDIF
    CASE('V')
       IF ( Vars_Existence( 'glamv' , parent_coordinate_file ) ) THEN
          lonParent = G0%glamv
          latParent = G0%gphiv
       ELSE
          WRITE(*,*) 'ERROR ***** missing glamv/gphiv coordinates in ', TRIM(parent_coordinate_file)
          STOP
       ENDIF
       IF ( Vars_Existence( 'glamv' , Childcoordinates ) ) THEN
          lonChild  = G1%glamv
          latChild  = G1%gphiv
       ELSE
          WRITE(*,*) 'ERROR ***** missing glamv/gphiv coordinates in ', TRIM(Childcoordinates)
          STOP
       ENDIF
       !!IF ( Vars_Existence( 'vmask' , Childcoordinates ) ) THEN
          mask = G1%vmask
       !!ELSE
       !!   mask = 1.
       !!ENDIF
    END SELECT

    DEALLOCATE(G0%glamu,G0%glamv,G0%gphiu,G0%gphiv)
    DEALLOCATE(G1%glamu,G1%glamv,G1%gphiu,G1%gphiv) 
    DEALLOCATE(G1%glamt,G1%gphit,G0%glamt,G0%gphit)

    !
    !longitude modification if child domain covers Pacific ocean area
    ! 
    IF( lonChild(1,1) > lonChild(nxfin,nyfin) ) THEN
       Pacifique = .TRUE.
       WHERE( lonChild < 0 )
          lonChild = lonChild + 360.
       END WHERE
       WHERE( lonParent < 0 )
          lonParent = lonParent + 360.
       END WHERE
    ENDIF

    !    
    !
    ! dimensions initialization      
    !
    CALL Write_Ncdf_dim('x',Child_file,nxfin)
    CALL Write_Ncdf_dim('y',Child_file,nyfin)
    IF ( Dims_Existence( 'deptht' , filename ) ) CALL Write_Ncdf_dim('deptht',Child_file,deptht)
    IF ( Dims_Existence( 'depthu' , filename ) ) CALL Write_Ncdf_dim('depthu',Child_file,deptht)
    IF ( Dims_Existence( 'depthv' , filename ) ) CALL Write_Ncdf_dim('depthv',Child_file,deptht)
    IF ( Dims_Existence( 'depthw' , filename ) ) CALL Write_Ncdf_dim('depthw',Child_file,deptht)
    IF ( Dims_Existence( 'z'      , filename ) ) CALL Write_Ncdf_dim('z'     ,Child_file,deptht)
    IF ( Dims_Existence( 'time_counter' , filename ) ) CALL Write_Ncdf_dim('time_counter',Child_file,time)
    IF ( Dims_Existence( 'time'         , filename ) ) CALL Write_Ncdf_dim('time'        ,Child_file,time)

    IF( deptht .NE. 1 .AND. deptht .NE. N ) THEN
       WRITE(*,*) '***'
       WRITE(*,*) 'Number of vertical levels doesn t match between namelist'
       WRITE(*,*) 'and forcing file ',TRIM(filename)
       WRITE(*,*) 'Please check the values in namelist file'
       WRITE(*,*) 'N = ',N
       WRITE(*,*) 'deptht = ',deptht
       STOP     
    ENDIF
    !
    !
    DO i = 1,SIZE(Ncdf_varname)
       !      
       !
       ! ******************************LOOP ON VARIABLE NAMES*******************************************
       !
       !      
       SELECT CASE (TRIM(Ncdf_varname(i)))
          !
          !copy nav_lon from child coordinates to output file      
       CASE('nav_lon')   
          ALLOCATE(latlon_temp(nxfin,nyfin))
          CALL Read_Ncdf_var('nav_lon',TRIM(Childcoordinates),latlon_temp) 
          CALL Write_Ncdf_var('nav_lon',(/'x','y'/),Child_file,latlon_temp,'float')
          CALL Copy_Ncdf_att('nav_lon',TRIM(filename),Child_file, &
               MINVAL(latlon_temp),MAXVAL(latlon_temp))
          DEALLOCATE(latlon_temp)
          varname = TRIM(Ncdf_varname(i))
          Interpolation = .FALSE.
          !	     
          !copy nav_lat from child coordinates to output file
       CASE('nav_lat')             
          ALLOCATE(latlon_temp(nxfin,nyfin))
          CALL Read_Ncdf_var('nav_lat',TRIM(Childcoordinates),latlon_temp) 
          CALL Write_Ncdf_var('nav_lat',(/'x','y'/),Child_file,latlon_temp,'float')
          CALL Copy_Ncdf_att('nav_lat',TRIM(filename),Child_file, &
               MINVAL(latlon_temp),MAXVAL(latlon_temp)) 
          DEALLOCATE(latlon_temp)
          varname = TRIM(Ncdf_varname(i))
          Interpolation = .FALSE.
          !
          !copy nav_lev from restart_file to output file
          !
       CASE('nav_lev')

          WRITE(*,*) 'copy nav_lev'
          CALL Read_Ncdf_var('nav_lev',filename,nav_lev) 
          IF(.NOT. dimg ) THEN
             CALL Write_Ncdf_var('nav_lev','z',Child_file,nav_lev,'float')
             CALL Copy_Ncdf_att('nav_lev',filename,Child_file)      
          ENDIF
          DEALLOCATE(nav_lev)
          Interpolation = .FALSE.
          !	    	
          !copy time_counter from input file to output file
       CASE('time_counter')
          ALLOCATE(timedepth_temp(time))
          CALL Read_Ncdf_var('time_counter',filename,timedepth_temp) 
          CALL Write_Ncdf_var('time_counter','time_counter',  &
               Child_file,timedepth_temp,'float')
          CALL Copy_Ncdf_att('time_counter',TRIM(filename),Child_file)
          DEALLOCATE(timedepth_temp)
          varname = TRIM(Ncdf_varname(i))
          Interpolation = .FALSE.
          !
          !copy time from input file to output file
       CASE('time')
          ALLOCATE(timedepth_temp(time))
          CALL Read_Ncdf_var('time',filename,timedepth_temp) 
          CALL Write_Ncdf_var('time','time',  &
               Child_file,timedepth_temp,'float')
          CALL Copy_Ncdf_att('time',TRIM(filename),Child_file)
          DEALLOCATE(timedepth_temp)
          varname = TRIM(Ncdf_varname(i))
          Interpolation = .FALSE.
          !
          !copy deptht from input file to output file
       CASE('deptht')
          ALLOCATE(depth(deptht))
          CALL Read_Ncdf_var('deptht',filename,depth) 
          CALL Write_Ncdf_var('deptht','deptht',Child_file,depth,'float')
          CALL Copy_Ncdf_att('deptht',TRIM(filename),Child_file)
          varname = TRIM(Ncdf_varname(i))
          Interpolation = .FALSE.     
          !
          !copy depthu from input file to output file
       CASE('depthu')
          ALLOCATE(depth(deptht))
          CALL Read_Ncdf_var('depthu',filename,depth) 
          CALL Write_Ncdf_var('depthu','depthu',Child_file,depth,'float')
          CALL Copy_Ncdf_att('depthu',TRIM(filename),Child_file)
          varname = TRIM(Ncdf_varname(i))
          Interpolation = .FALSE.     
          !
          !copy depthv from input file to output file
       CASE('depthv')
          ALLOCATE(depth(deptht))
          CALL Read_Ncdf_var('depthv',filename,depth) 
          CALL Write_Ncdf_var('depthv','depthv',Child_file,depth,'float')
          CALL Copy_Ncdf_att('depthv',TRIM(filename),Child_file)
          varname = TRIM(Ncdf_varname(i))
          Interpolation = .FALSE.     
          !
          !copy depthv from input file to output file
       CASE('depthw')
          ALLOCATE(depth(deptht))
          CALL Read_Ncdf_var('depthw',filename,depth) 
          CALL Write_Ncdf_var('depthw','depthw',Child_file,depth,'float')
          CALL Copy_Ncdf_att('depthw',TRIM(filename),Child_file)
          varname = TRIM(Ncdf_varname(i))
          Interpolation = .FALSE.     
          !
          !copy time_steps from input file in case of restget use in NEMO in/out routines
       CASE('time_steps')
          !         print *,'avant time step'

          CALL Read_Ncdf_var('time_steps',filename,tabtemp1DInt) 
          !		   print *,'timedeph = ',tabtemp1DInt
          CALL Write_Ncdf_var('time_steps','time_counter',Child_file,tabtemp1DInt,'integer')
          CALL Copy_Ncdf_att('time_steps',filename,Child_file) 
          DEALLOCATE(tabtemp1DInt)
          Interpolation = .FALSE.  
          !  
          !store tmask in output file
       CASE('tmask')
          dimnames(1)='x'
          dimnames(2)='y'
          dimnames(3)='deptht' 
          IF (.NOT.ASSOCIATED(G1%tmask)) THEN
             ALLOCATE(G1%tmask(nxfin,nyfin,deptht))
             G1%tmask = 1
          ENDIF
          CALL Write_Ncdf_var('tmask',dimnames(1:3),Child_file,G1%tmask,'float')
          varname = TRIM(Ncdf_varname(i))
          Interpolation = .FALSE.     
          !      
       CASE default
          varname = Ncdf_varname(i)
          conservation = .FALSE.
          CALL get_interptype( varname,interp_type,conservation )
          WRITE(*,*) '**********************************************'
          WRITE(*,*) 'varname      = ',TRIM(varname), 'at ', posvar, ' point'
          WRITE(*,*) 'interp_type  = ',TRIM(interp_type)
          WRITE(*,*) 'conservation = ',conservation
          WRITE(*,*) '***********************************************'
          Interpolation = .TRUE.
          !	         
       END SELECT

       ! //////////////// INTERPOLATION FOR 2D VARIABLES /////////////////////////////////////
       ! 
       IF( Interpolation .AND. Get_NbDims(TRIM(varname),TRIM(filename)) == 2 ) THEN
          !      
          ALLOCATE(detected_pts(numlon,numlat,N))
          ALLOCATE(masksrc(numlon,numlat))                                           
          !
          IF(extrapolation) THEN
             WRITE(*,*) 'interpolation/extrapolation ',TRIM(varname)
          ELSE
             WRITE(*,*) 'interpolation ',TRIM(varname)
          ENDIF
          !                            
          ALLOCATE(tabvar3d(numlon,numlat,1))      
          ALLOCATE(tabinterp3d(nxfin,nyfin,1))
          !                    
          CALL Read_Ncdf_var(varname,filename,tabvar3d)                                              
          !
          ! search points where extrapolation is required
          ! 
          IF(Extrapolation) THEN
             WHERE( tabvar3d .GE. 1.e+20 ) tabvar3d = 0.
             CALL extrap_detect(G0,G1,detected_pts(:,:,1),1)
             CALL correct_field_2d(detected_pts(:,:,1),tabvar3d,G0,masksrc,'posvar')
          ELSE
             masksrc = .TRUE.
          ENDIF

          SELECT CASE(TRIM(interp_type))
          CASE('bilinear')
             CALL get_remap_matrix(latParent,latChild,   &
                lonParent,lonChild,                    &
                masksrc,matrix,src_add,dst_add)

          CASE('bicubic')
             CALL get_remap_bicub(latParent,latChild,   &
                lonParent,lonChild,   &
                masksrc,matrix,src_add,dst_add)

          END SELECT
          !                                                        
          !      
          SELECT CASE(TRIM(interp_type))
          CASE('bilinear')                                                       
             CALL make_remap(tabvar3d(:,:,1),tabinterp3d(:,:,1),nxfin,nyfin, &
                matrix,src_add,dst_add)     
          CASE('bicubic')                                   
             CALL make_bicubic_remap(tabvar3d(:,:,1),masksrc,tabinterp3d(:,:,1),nxfin,nyfin, &
                matrix,src_add,dst_add)                        
          END SELECT
          !                     
          IF( conservation ) CALL Correctforconservation(tabvar3d(:,:,1),tabinterp3d(:,:,1), &
             G0%e1t,G0%e2t,G1%e1t,G1%e2t,nxfin,nyfin,posvar,imin,jmin)  
          !      
          IF(Extrapolation) tabinterp3d(:,:,1) = tabinterp3d(:,:,1) * mask(:,:,1)     
!!          tabinterp3d(:,:,1) = tabinterp3d(:,:,1) * mask(:,:,1)
          !                     
          dimnames(1)='x'
          dimnames(2)='y'
          !                                     
          CALL Write_Ncdf_var(TRIM(varname),dimnames(1:2),&
             Child_file,tabinterp3d(:,:,1),'float')
          !      
          DEALLOCATE(tabinterp3d)
          DEALLOCATE(tabvar3d)                     
          !                   
          DEALLOCATE(detected_pts)
          IF(ASSOCIATED(matrix)) DEALLOCATE(matrix,dst_add,src_add)
          DEALLOCATE( masksrc)
          
          CALL Copy_Ncdf_att(TRIM(varname),TRIM(filename),Child_file)                        
          !
          
       ! //////////////// INTERPOLATION FOR 3D VARIABLES /////////////////////////////////////
       ! 
       ELSEIF( Interpolation .AND. Get_NbDims(TRIM(varname),TRIM(filename)) == 3 ) THEN
          !      
          IF( DimUnlimited_Var(TRIM(varname),TRIM(filename)) ) THEN

             ALLOCATE(detected_pts(numlon,numlat,N))
             ALLOCATE(masksrc(numlon,numlat))                                           
             !
             ! ******************************LOOP ON TIME*******************************************
             !loop on time
             DO t = 1,time
                !                   
                IF(extrapolation) THEN
                   WRITE(*,*) 'interpolation/extrapolation ',TRIM(varname),' for time t = ',t
                ELSE
                   WRITE(*,*) 'interpolation ',TRIM(varname),' for time t = ',t
                ENDIF
                !                            
                ALLOCATE(tabvar3d(numlon,numlat,1))      
                ALLOCATE(tabinterp3d(nxfin,nyfin,1))
                !                    
                CALL Read_Ncdf_var(varname,filename,tabvar3d,t)                                              
                !
                ! search points where extrapolation is required
                ! 
                IF(Extrapolation) THEN
                   WHERE( tabvar3d .GE. 1.e+20 ) tabvar3d = 0.
                   IF (t .EQ. 1. ) CALL extrap_detect(G0,G1,detected_pts(:,:,1),1)
                   CALL correct_field_2d(detected_pts(:,:,1),tabvar3d,G0,masksrc,'posvar')
                ELSE
                   masksrc = .TRUE.
                ENDIF

                IF (t.EQ.1 ) THEN 

                   SELECT CASE(TRIM(interp_type))
                   CASE('bilinear')
                      CALL get_remap_matrix(latParent,latChild,   &
                         lonParent,lonChild,                    &
                         masksrc,matrix,src_add,dst_add)
                   CASE('bicubic')
                      CALL get_remap_bicub(latParent,latChild,   &
                         lonParent,lonChild,   &
                         masksrc,matrix,src_add,dst_add)

                   END SELECT
                   !                                                        
                ENDIF
                !      
                SELECT CASE(TRIM(interp_type))
                CASE('bilinear')                                                       
                   CALL make_remap(tabvar3d(:,:,1),tabinterp3d(:,:,1),nxfin,nyfin, &
                      matrix,src_add,dst_add)     
                CASE('bicubic')                                   
                   CALL make_bicubic_remap(tabvar3d(:,:,1),masksrc,tabinterp3d(:,:,1),nxfin,nyfin, &
                      matrix,src_add,dst_add)                        
                END SELECT
                !                     
                IF( conservation ) CALL Correctforconservation(tabvar3d(:,:,1),tabinterp3d(:,:,1), &
                   G0%e1t,G0%e2t,G1%e1t,G1%e2t,nxfin,nyfin,posvar,imin,jmin)  
                !      
                IF(Extrapolation) tabinterp3d(:,:,1) = tabinterp3d(:,:,1) * mask(:,:,1)     
!!                tabinterp3d(:,:,1) = tabinterp3d(:,:,1) * mask(:,:,1)
                !                     
                dimnames(1)='x'
                dimnames(2)='y'
                dimnames(3)='time_counter'
                !                                     
                CALL Write_Ncdf_var(TRIM(varname),dimnames(1:3),&
                   Child_file,tabinterp3d,t,'float')
                !      
                DEALLOCATE(tabinterp3d)
                DEALLOCATE(tabvar3d)                     
                !end loop on time
             END DO
             !                   
             DEALLOCATE(detected_pts)
             IF(ASSOCIATED(matrix)) DEALLOCATE(matrix,dst_add,src_add)
             DEALLOCATE( masksrc)

          ELSE

             dimnames(1)='x'
             dimnames(2)='y'
             IF ( Dims_Existence( 'deptht' , filename ) ) dimnames(3)='deptht'  
             IF ( Dims_Existence( 'depthu' , filename ) ) dimnames(3)='depthu'
             IF ( Dims_Existence( 'depthv' , filename ) ) dimnames(3)='depthv'
             IF ( Dims_Existence( 'depthw' , filename ) ) dimnames(3)='depthw'
             IF ( Dims_Existence( 'z' , filename ) ) dimnames(3)='z'
             !
             ! loop on vertical levels
             ! 
             DO nb = 1,deptht
                !        
                ALLOCATE(masksrc(numlon,numlat))
                ALLOCATE(detected_pts(numlon,numlat,N))          
                !
                ! point detection et level n                                     
                !
                land_level = .FALSE.
                IF( Extrapolation ) THEN
                   IF(MAXVAL(mask(:,:,nb))==0.) land_level = .TRUE.
                ENDIF


                IF ( land_level ) THEN
                   !
                   WRITE(*,*) 'only land points on level ',nb
                   ALLOCATE(tabinterp3d(nxfin,nyfin,1)) 
                   tabinterp3d = 0.e0
                   !                              
                   CALL Write_Ncdf_var(TRIM(varname),dimnames, &
                      Child_file,tabinterp3d,nb,'float')
                   DEALLOCATE(tabinterp3d)
                   !
                ELSE
                   !                       
                   ALLOCATE(tabvar01(numlon,numlat,1))    ! level k
                   IF(Extrapolation) ALLOCATE(tabvar02(numlon,numlat,1))    ! level k-1
                   IF(Extrapolation) ALLOCATE(tabvar03(numlon,numlat,1))    ! level k-2                       
                   ALLOCATE(tabinterp3d(nxfin,nyfin,1)) 
                   !
                   IF(Extrapolation) THEN                                  
                      IF(nb==1) THEN
                         CALL Read_Ncdf_var(varname,filename,tabvar01,nb)    
                         WHERE( tabvar01 .GE. 1.e+20 ) tabvar01 = 0.
                      ELSE IF (nb==2) THEN
                         CALL Read_Ncdf_var(varname,filename,tabvar02,nb-1)
                         CALL Read_Ncdf_var(varname,filename,tabvar01,nb)            
                         WHERE( tabvar01 .GE. 1.e+20 ) tabvar01 = 0.
                         WHERE( tabvar02 .GE. 1.e+20 ) tabvar02 = 0.
                      ELSE 
                         CALL Read_Ncdf_var(varname,filename,tabvar03,nb-2)
                         CALL Read_Ncdf_var(varname,filename,tabvar02,nb-1)
                         CALL Read_Ncdf_var(varname,filename,tabvar01,nb)
                         WHERE( tabvar01 .GE. 1.e+20 ) tabvar01 = 0.
                         WHERE( tabvar02 .GE. 1.e+20 ) tabvar02 = 0.
                         WHERE( tabvar03 .GE. 1.e+20 ) tabvar03 = 0.
                      ENDIF
                      !                              
                      CALL extrap_detect(G0,G1,detected_pts(:,:,nb),nb)

                      ALLOCATE(tabvar1(numlon,numlat,1,1))    ! level k
                      ALLOCATE(tabvar2(numlon,numlat,1,1))    ! level k-1
                      ALLOCATE(tabvar3(numlon,numlat,1,1))    ! level k-2
                      tabvar1(:,:,1,1) = tabvar01(:,:,1)
                      tabvar2(:,:,1,1) = tabvar02(:,:,1)
                      tabvar3(:,:,1,1) = tabvar03(:,:,1)
                      CALL correct_field(detected_pts(:,:,nb),tabvar1,tabvar2,&
                         tabvar3,G0,depth,masksrc,nb,'posvar')
                      tabvar01(:,:,1) = tabvar1(:,:,1,1)
                      tabvar02(:,:,1) = tabvar2(:,:,1,1)
                      tabvar03(:,:,1) = tabvar3(:,:,1,1)                      
                      DEALLOCATE(tabvar1,tabvar2,tabvar3)
                      DEALLOCATE(tabvar02,tabvar03)

                   ELSE 
                      CALL Read_Ncdf_var(varname,filename,tabvar01,nb)
                      IF(MAXVAL(tabvar01(:,:,1))==0.) land_level = .TRUE.
                      masksrc = .TRUE. 
                   ENDIF

                   IF( Extrapolation ) THEN
                      WRITE(*,*) 'interpolation/extrapolation ',TRIM(varname),' for vertical level = ',nb   
                   ELSE
                      WRITE(*,*) 'interpolation ',TRIM(varname),' for vertical level = ',nb
                   ENDIF
                   !
                   SELECT CASE(TRIM(interp_type))
                   CASE('bilinear')
                      CALL get_remap_matrix(latParent,latChild,   &
                         lonParent,lonChild,   &
                         masksrc,matrix,src_add,dst_add)

                   CASE('bicubic')
                      CALL get_remap_bicub(latParent,latChild,   &
                         lonParent,lonChild,   &
                         masksrc,matrix,src_add,dst_add)
                      !                             
                   END SELECT
                   !                                                        
                   !      
                   SELECT CASE(TRIM(interp_type))
                      !                             
                   CASE('bilinear')                                                       
                      CALL make_remap(tabvar01(:,:,1),tabinterp3d(:,:,1),nxfin,nyfin, &
                         matrix,src_add,dst_add)     
                   CASE('bicubic')                                   
                      CALL make_bicubic_remap(tabvar01(:,:,1),masksrc,tabinterp3d(:,:,1),nxfin,nyfin, &
                         matrix,src_add,dst_add)                        
                   END SELECT
                   !                     
                   IF( conservation ) CALL Correctforconservation(tabvar01(:,:,1),tabinterp3d(:,:,1), &
                      G0%e1t,G0%e2t,G1%e1t,G1%e2t,nxfin,nyfin,posvar,imin,jmin)  

                   !
                   ALLOCATE(tabinterp4d(nxfin,nyfin,1,1))
                   tabinterp4d(:,:,1,1) = tabinterp3d(:,:,1)
                   IF(Extrapolation) CALL check_extrap(G1,tabinterp4d,nb)
                   tabinterp3d(:,:,1) = tabinterp4d(:,:,1,1)
                   DEALLOCATE(tabinterp4d)
                   !     
                   IF(Extrapolation) tabinterp3d(:,:,1) = tabinterp3d(:,:,1) * mask(:,:,nb)     
!!                   tabinterp3d(:,:,1) = tabinterp3d(:,:,1) * mask(:,:,nb)
                   !                     
                   CALL Write_Ncdf_var(TRIM(varname),dimnames, &
                      Child_file,tabinterp3d,nb,'float')
                   !      
                   DEALLOCATE(tabinterp3d)
                   DEALLOCATE(tabvar01)                     
                   !
                   !
                ENDIF

                !
                IF(ASSOCIATED(matrix)) DEALLOCATE(matrix,dst_add,src_add)
                DEALLOCATE( masksrc )
                DEALLOCATE(detected_pts)       
                !
                ! end loop on vertical levels
                !      
             END DO

             
          ENDIF

          CALL Copy_Ncdf_att(TRIM(varname),TRIM(filename),Child_file)                        
          !
       ELSE IF( Interpolation .AND. Get_NbDims(TRIM(varname),TRIM(filename)) == 4 ) THEN
          !
          !
          ! //////////////// INTERPOLATION FOR 4D VARIABLES /////////////////////////////////////
          ! 
          dimnames(1)='x'
          dimnames(2)='y'
          IF ( Dims_Existence( 'deptht' , filename ) ) dimnames(3)='deptht'  
          IF ( Dims_Existence( 'depthu' , filename ) ) dimnames(3)='depthu'
          IF ( Dims_Existence( 'depthv' , filename ) ) dimnames(3)='depthv'
          IF ( Dims_Existence( 'depthw' , filename ) ) dimnames(3)='depthw'
          IF ( Dims_Existence( 'z' , filename ) ) dimnames(3)='z'
          dimnames(4)='time_counter'  
          !
          ! loop on vertical levels
          ! 
          DO nb = 1,deptht
             !        
             ALLOCATE(masksrc(numlon,numlat))
             ALLOCATE(detected_pts(numlon,numlat,N))          
             !
             ! point detection et level n                                     
             !
             land_level = .FALSE.
             IF( Extrapolation ) THEN
                IF(MAXVAL(mask(:,:,nb))==0.) land_level = .TRUE.
             ENDIF


             IF ( land_level ) THEN
                !
                WRITE(*,*) 'only land points on level ',nb
                ALLOCATE(tabinterp4d(nxfin,nyfin,1,1)) 
                tabinterp4d = 0.e0
                !                              
                DO ii = 1,time
                   CALL Write_Ncdf_var(TRIM(varname),dimnames, &
                        Child_file,tabinterp4d,ii,nb,'float')
                END DO
                DEALLOCATE(tabinterp4d)
                !
             ELSE
                !                       
                ! loop on time                       
                !
                DO t = 1,time
                   !                   
                   ALLOCATE(tabvar1(numlon,numlat,1,1))    ! level k
                   IF(Extrapolation) ALLOCATE(tabvar2(numlon,numlat,1,1))    ! level k-1
                   IF(Extrapolation) ALLOCATE(tabvar3(numlon,numlat,1,1))    ! level k-2                       
                   ALLOCATE(tabinterp4d(nxfin,nyfin,1,1)) 
                   !
                   IF(Extrapolation) THEN                                  
                      IF(nb==1) THEN
                         CALL Read_Ncdf_var(varname,filename,tabvar1,t,nb)    
                         WHERE( tabvar1 .GE. 1.e+20 ) tabvar1 = 0.
                      ELSE IF (nb==2) THEN
                         CALL Read_Ncdf_var(varname,filename,tabvar2,t,nb-1)
                         CALL Read_Ncdf_var(varname,filename,tabvar1,t,nb)            
                         WHERE( tabvar1 .GE. 1.e+20 ) tabvar1 = 0.
                         WHERE( tabvar2 .GE. 1.e+20 ) tabvar2 = 0.
                      ELSE 
                         CALL Read_Ncdf_var(varname,filename,tabvar3,t,nb-2)
                         CALL Read_Ncdf_var(varname,filename,tabvar2,t,nb-1)
                         CALL Read_Ncdf_var(varname,filename,tabvar1,t,nb)
                         WHERE( tabvar1 .GE. 1.e+20 ) tabvar1 = 0.
                         WHERE( tabvar2 .GE. 1.e+20 ) tabvar2 = 0.
                         WHERE( tabvar3 .GE. 1.e+20 ) tabvar3 = 0.
                      ENDIF
                      !                              
                      IF (t.EQ.1 ) CALL extrap_detect(G0,G1,detected_pts(:,:,nb),nb)

                      CALL correct_field(detected_pts(:,:,nb),tabvar1,tabvar2,&
                           tabvar3,G0,depth,masksrc,nb,'posvar')                       
                      DEALLOCATE(tabvar2,tabvar3)

                   ELSE 
                      CALL Read_Ncdf_var(varname,filename,tabvar1,t,nb)
                      IF(MAXVAL(tabvar1(:,:,1,1))==0.) land_level = .TRUE.
                      masksrc = .TRUE. 
                   ENDIF

                   IF( Extrapolation ) THEN
                      WRITE(*,*) 'interpolation/extrapolation ',TRIM(varname),' for time t = ',t,'vertical level = ',nb   
                   ELSE
                      WRITE(*,*) 'interpolation ',TRIM(varname),' for time t = ',t,'vertical level = ',nb
                   ENDIF

                   !                          

                   IF (t.EQ.1 ) THEN 

                      SELECT CASE(TRIM(interp_type))
                      CASE('bilinear')
                         CALL get_remap_matrix(latParent,latChild,   &
                              lonParent,lonChild,   &
                              masksrc,matrix,src_add,dst_add)

                      CASE('bicubic')
                         CALL get_remap_bicub(latParent,latChild,   &
                              lonParent,lonChild,   &
                              masksrc,matrix,src_add,dst_add)
                         !                             
                      END SELECT
                      !                                                        
                   ENDIF
                   !      
                   SELECT CASE(TRIM(interp_type))
                      !                             
                   CASE('bilinear')                                                       
                      CALL make_remap(tabvar1(:,:,1,1),tabinterp4d(:,:,1,1),nxfin,nyfin, &
                           matrix,src_add,dst_add)     
                   CASE('bicubic')                                   
                      CALL make_bicubic_remap(tabvar1(:,:,1,1),masksrc,tabinterp4d(:,:,1,1),nxfin,nyfin, &
                           matrix,src_add,dst_add)                        
                   END SELECT
                   !                     
                   IF( conservation ) CALL Correctforconservation(tabvar1(:,:,1,1),tabinterp4d(:,:,1,1), &
                        G0%e1t,G0%e2t,G1%e1t,G1%e2t,nxfin,nyfin,posvar,imin,jmin)  

                   !
                   IF(Extrapolation) CALL check_extrap(G1,tabinterp4d,nb)
                   !      
                   IF(Extrapolation) tabinterp4d(:,:,1,1) = tabinterp4d(:,:,1,1) * mask(:,:,nb)     
                   !                     
                   CALL Write_Ncdf_var(TRIM(varname),dimnames, &
                        Child_file,tabinterp4d,t,nb,'float')
                   !      
                   DEALLOCATE(tabinterp4d)
                   DEALLOCATE(tabvar1)                     
                   !
                   ! end loop on time
                   !

                END DO

             ENDIF

             !
             IF(ASSOCIATED(matrix)) DEALLOCATE(matrix,dst_add,src_add)
             DEALLOCATE( masksrc )
             DEALLOCATE(detected_pts)       
             !
             ! end loop on vertical levels
             !      
          END DO
          !   
          CALL Copy_Ncdf_att(TRIM(varname),TRIM(filename),Child_file)
          !
          ! fin du if interpolation ...
          !
       ENDIF
       ! 
    END DO

    PRINT *,'FIN DE INTERPEXTRAPVAR'
    !    
    IF(Extrapolation) DEALLOCATE(G1%tmask,G0%tmask)
    DEALLOCATE(G0%e1t,G0%e2t,G1%e1t,G1%e2t)
    IF(ASSOCIATED(depth)) DEALLOCATE(depth)   
    !
  END SUBROUTINE Interp_Extrap_var
  !
  ! 
END MODULE agrif_interpolation
