












      SUBROUTINE linterp2d (LBx, UBx, LBy, UBy,
     &                      Xinp, Yinp, Finp,
     &                      Istr, Iend, Jstr, Jend,
     &                      Xout, Yout, Fout)      


      use Agrif_Util
      interface
        subroutine Sub_Loop_linterp2d(LBx,UBx,LBy,UBy,Xinp,Yinp,Finp,Ist
     &r,Iend,Jstr,Jend,Xout,Yout,Fout,padd_E,Mm,padd_X,Lm)
          implicit none
      integer(4), parameter :: LLm0 = 41
      integer(4), parameter :: MMm0 = 42
      integer(4), parameter :: N = 32
      integer(4), parameter :: NP_XI = 1
      integer(4), parameter :: NP_ETA = 4
      integer(4), parameter :: NNODES = NP_XI*NP_ETA
      integer(4), parameter :: NPP = 1
      integer(4), parameter :: NSUB_X = 1
      integer(4), parameter :: NSUB_E = 1
      integer(4), parameter :: NWEIGHT = 1000
      integer(4), parameter :: stdout = 6
      integer(4), parameter :: Np = N+1
      integer(4), parameter :: NSA = 28
      real, parameter :: Vtransform = 2
      integer(4), parameter :: itemp = 1
      integer(4), parameter :: ntrc_temp = 1
      integer(4), parameter :: ntrc_salt = 1
      integer(4), parameter :: ntrc_pas = 0
      integer(4), parameter :: ntrc_bio = 0
      integer(4), parameter :: ntrc_subs = 0
      integer(4), parameter :: ntrc_substot = 0
      integer(4), parameter :: ntrc_sed = 0
      integer(4), parameter :: NTA = itemp+ntrc_salt
      integer(4), parameter :: NT = itemp+ntrc_salt+ntrc_pas+ntrc_bio+nt
     &rc_sed
      integer(4), parameter :: NTot = NT
      integer(4), parameter :: isalt = itemp+1
      integer(4), parameter :: ntrc_diabio = 0
      integer(4), parameter :: ntrc_diats = 0
      integer(4), parameter :: ntrc_diauv = 0
      integer(4), parameter :: ntrc_diavrt = 0
      integer(4), parameter :: ntrc_diaek = 0
      integer(4), parameter :: ntrc_diapv = 0
      integer(4), parameter :: ntrc_diaeddy = 0
      integer(4), parameter :: ntrc_surf = 0
      integer(4) :: padd_E
      integer(4) :: Mm
      integer(4) :: padd_X
      integer(4) :: Lm
      integer(4), intent(in) :: LBx
      integer(4), intent(in) :: UBx
      integer(4), intent(in) :: LBy
      integer(4), intent(in) :: UBy
      integer(4), intent(in) :: Istr
      integer(4), intent(in) :: Iend
      integer(4), intent(in) :: Jstr
      integer(4), intent(in) :: Jend
      real(8), intent(in), dimension(LBx:UBx) :: Xinp
      real(8), intent(in), dimension(LBy:UBy) :: Yinp
      real(8), intent(in), dimension(LBx:UBx,LBy:UBy) :: Finp
      real(8), intent(in), dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: X
     &out
      real(8), intent(in), dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: Y
     &out
      real(8), intent(out), dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: 
     &Fout
        end subroutine Sub_Loop_linterp2d

      end interface
      integer(4), intent(in) :: LBx
      integer(4), intent(in) :: UBx
      integer(4), intent(in) :: LBy
      integer(4), intent(in) :: UBy
      integer(4), intent(in) :: Istr
      integer(4), intent(in) :: Iend
      integer(4), intent(in) :: Jstr
      integer(4), intent(in) :: Jend
      real(8), intent(in), dimension(LBx:UBx) :: Xinp
      real(8), intent(in), dimension(LBy:UBy) :: Yinp
      real(8), intent(in), dimension(LBx:UBx,LBy:UBy) :: Finp
      real(8), intent(in), dimension(-1: Agrif_tabvars_i(2)%iarray0+2+ A
     &grif_tabvars_i(13)%iarray0,-1: Agrif_tabvars_i(4)%iarray0+2+ Agrif
     &_tabvars_i(14)%iarray0) :: Xout
      real(8), intent(in), dimension(-1: Agrif_tabvars_i(2)%iarray0+2+ A
     &grif_tabvars_i(13)%iarray0,-1: Agrif_tabvars_i(4)%iarray0+2+ Agrif
     &_tabvars_i(14)%iarray0) :: Yout
      real(8), intent(out), dimension(-1: Agrif_tabvars_i(2)%iarray0+2+ 
     &Agrif_tabvars_i(13)%iarray0,-1: Agrif_tabvars_i(4)%iarray0+2+ Agri
     &f_tabvars_i(14)%iarray0) :: Fout
      

        call Sub_Loop_linterp2d(LBx,UBx,LBy,UBy,Xinp,Yinp,Finp,Istr,Iend
     &,Jstr,Jend,Xout,Yout,Fout, Agrif_tabvars_i(14)%iarray0, Agrif_tabv
     &ars_i(4)%iarray0, Agrif_tabvars_i(13)%iarray0, Agrif_tabvars_i(2)%
     &iarray0)

      end


      subroutine Sub_Loop_linterp2d(LBx,UBx,LBy,UBy,Xinp,Yinp,Finp,Istr,
     &Iend,Jstr,Jend,Xout,Yout,Fout,padd_E,Mm,padd_X,Lm)

      use Agrif_Types, only : Agrif_tabvars

      implicit none

      integer(4), parameter :: LLm0 = 41
      integer(4), parameter :: MMm0 = 42
      integer(4), parameter :: N = 32
      integer(4), parameter :: NP_XI = 1
      integer(4), parameter :: NP_ETA = 4
      integer(4), parameter :: NNODES = NP_XI*NP_ETA
      integer(4), parameter :: NPP = 1
      integer(4), parameter :: NSUB_X = 1
      integer(4), parameter :: NSUB_E = 1
      integer(4), parameter :: NWEIGHT = 1000
      integer(4), parameter :: stdout = 6
      integer(4), parameter :: Np = N+1
      integer(4), parameter :: NSA = 28
      real, parameter :: Vtransform = 2
      integer(4), parameter :: itemp = 1
      integer(4), parameter :: ntrc_temp = 1
      integer(4), parameter :: ntrc_salt = 1
      integer(4), parameter :: ntrc_pas = 0
      integer(4), parameter :: ntrc_bio = 0
      integer(4), parameter :: ntrc_subs = 0
      integer(4), parameter :: ntrc_substot = 0
      integer(4), parameter :: ntrc_sed = 0
      integer(4), parameter :: NTA = itemp+ntrc_salt
      integer(4), parameter :: NT = itemp+ntrc_salt+ntrc_pas+ntrc_bio+nt
     &rc_sed
      integer(4), parameter :: NTot = NT
      integer(4), parameter :: isalt = itemp+1
      integer(4), parameter :: ntrc_diabio = 0
      integer(4), parameter :: ntrc_diats = 0
      integer(4), parameter :: ntrc_diauv = 0
      integer(4), parameter :: ntrc_diavrt = 0
      integer(4), parameter :: ntrc_diaek = 0
      integer(4), parameter :: ntrc_diapv = 0
      integer(4), parameter :: ntrc_diaeddy = 0
      integer(4), parameter :: ntrc_surf = 0
      integer(4) :: padd_E
      integer(4) :: Mm
      integer(4) :: padd_X
      integer(4) :: Lm
      integer(4), intent(in) :: LBx
      integer(4), intent(in) :: UBx
      integer(4), intent(in) :: LBy
      integer(4), intent(in) :: UBy
      integer(4), intent(in) :: Istr
      integer(4), intent(in) :: Iend
      integer(4), intent(in) :: Jstr
      integer(4), intent(in) :: Jend
      real(8), intent(in), dimension(LBx:UBx) :: Xinp
      real(8), intent(in), dimension(LBy:UBy) :: Yinp
      real(8), intent(in), dimension(LBx:UBx,LBy:UBy) :: Finp
      real(8), intent(in), dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: X
     &out
      real(8), intent(in), dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: Y
     &out
      real(8), intent(out), dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: 
     &Fout

                   

                                            
                             
                                            

                                                             
                                                            
                                           
                                                               
                                    
                                      
                                                          

                       

                                    

                        
                              

                                          
                                             
                          

                        

                                          
                        

                                                        
                      
                              

                                       
                                                                     
                                          
                         

                             

                             

                            

                            

                                             

                            

                                     

                                                               

                         

                                                      
                                                      
                                                                  
                               

                               

                              

                              

                               

                              

                              

                                

                             

                                                  
                                                      
                                                
                                                
                                                        
                                                                      
                                                                      
                                                                       
                                             
      integer(4) :: i
      integer(4) :: i1
      integer(4) :: i2
      integer(4) :: j
      integer(4) :: j1
      integer(4) :: j2
      integer(4) :: ii
      integer(4) :: jj
                                             
      real(8) :: cff
      real(8) :: x
      real(8) :: x1
      real(8) :: x2
      real(8) :: y
      real(8) :: y1
      real(8) :: y2
      DO j=Jstr,Jend
        DO i=Istr,Iend
           DO ii=LBx,(UBx-1)
             if ((Xinp(ii).le.Xout(i,j)).and.
     &           (Xinp(ii+1).gt.Xout(i,j))) then
               i1=ii
               i2=ii+1
               goto 10
             endif
           enddo
           print*, 'Did not find i1 and i2'
           goto 100
 10        continue
           DO jj=LBy,(UBy-1)
             if ((Yinp(jj).le.Yout(i,j)).and.
     &           (Yinp(jj+1).gt.Yout(i,j))) then
               j1=jj
               j2=jj+1
               goto 20
             endif
           enddo
           print*, 'Did not find j1 and j2'
           goto 100
 20        continue
          IF (((LBx.le.i1).and.(i1.le.UBx)).and.
     &        ((LBy.le.j1).and.(j1.le.UBy))) THEN
            x1=Xinp(i1)
            x2=Xinp(i2)
            y1=Yinp(j1)
            y2=Yinp(j2)
            x=Xout(i,j)
            y=Yout(i,j)
            cff= Finp(i1,j1)*(x2-x )*(y2-y )
     &          +Finp(i2,j1)*(x -x1)*(y2-y )
     &          +Finp(i1,j2)*(x2-x )*(y -y1)
     &          +Finp(i2,j2)*(x -x1)*(y -y1)
            Fout(i,j)=cff/((x2-x1)*(y2-y1))
          END IF
        END DO
      END DO
      RETURN
 100  continue
      print*, 'error in linterp2d'
             

      end subroutine Sub_Loop_linterp2d


      SUBROUTINE cinterp2d (LBx, UBx, LBy, UBy,
     &                      Xinp, Yinp, Finp,
     &                      Istr, Iend, Jstr, Jend,
     &                      Xout, Yout, Fout)      


      use Agrif_Util
      interface
        subroutine Sub_Loop_cinterp2d(LBx,UBx,LBy,UBy,Xinp,Yinp,Finp,Ist
     &r,Iend,Jstr,Jend,Xout,Yout,Fout,padd_E,Mm,padd_X,Lm)
          implicit none
      integer(4), parameter :: LLm0 = 41
      integer(4), parameter :: MMm0 = 42
      integer(4), parameter :: N = 32
      integer(4), parameter :: NP_XI = 1
      integer(4), parameter :: NP_ETA = 4
      integer(4), parameter :: NNODES = NP_XI*NP_ETA
      integer(4), parameter :: NPP = 1
      integer(4), parameter :: NSUB_X = 1
      integer(4), parameter :: NSUB_E = 1
      integer(4), parameter :: NWEIGHT = 1000
      integer(4), parameter :: stdout = 6
      integer(4), parameter :: Np = N+1
      integer(4), parameter :: NSA = 28
      real, parameter :: Vtransform = 2
      integer(4), parameter :: itemp = 1
      integer(4), parameter :: ntrc_temp = 1
      integer(4), parameter :: ntrc_salt = 1
      integer(4), parameter :: ntrc_pas = 0
      integer(4), parameter :: ntrc_bio = 0
      integer(4), parameter :: ntrc_subs = 0
      integer(4), parameter :: ntrc_substot = 0
      integer(4), parameter :: ntrc_sed = 0
      integer(4), parameter :: NTA = itemp+ntrc_salt
      integer(4), parameter :: NT = itemp+ntrc_salt+ntrc_pas+ntrc_bio+nt
     &rc_sed
      integer(4), parameter :: NTot = NT
      integer(4), parameter :: isalt = itemp+1
      integer(4), parameter :: ntrc_diabio = 0
      integer(4), parameter :: ntrc_diats = 0
      integer(4), parameter :: ntrc_diauv = 0
      integer(4), parameter :: ntrc_diavrt = 0
      integer(4), parameter :: ntrc_diaek = 0
      integer(4), parameter :: ntrc_diapv = 0
      integer(4), parameter :: ntrc_diaeddy = 0
      integer(4), parameter :: ntrc_surf = 0
      real(8), parameter :: C01 = 1.0D0/48.0D0
      real(8), parameter :: C02 = 1.0D0/32.0D0
      real(8), parameter :: C03 = 0.0625D0
      real(8), parameter :: C04 = 1.0D0/6.0D0
      real(8), parameter :: C05 = 0.25D0
      real(8), parameter :: C06 = 0.5D0
      real(8), parameter :: C07 = 0.3125D0
      real(8), parameter :: C08 = 0.625D0
      real(8), parameter :: C09 = 1.5D0
      real(8), parameter :: C10 = 13.0D0/24.0D0
      real(8), parameter :: LIMTR = 3.0D0
      real(8), parameter :: spv = 0.0D0
      integer(4) :: padd_E
      integer(4) :: Mm
      integer(4) :: padd_X
      integer(4) :: Lm
      integer(4), intent(in) :: LBx
      integer(4), intent(in) :: UBx
      integer(4), intent(in) :: LBy
      integer(4), intent(in) :: UBy
      integer(4), intent(in) :: Istr
      integer(4), intent(in) :: Iend
      integer(4), intent(in) :: Jstr
      integer(4), intent(in) :: Jend
      real(8), intent(in), dimension(LBx:UBx) :: Xinp
      real(8), intent(in), dimension(LBy:UBy) :: Yinp
      real(8), intent(in), dimension(LBx:UBx,LBy:UBy) :: Finp
      real(8), intent(in), dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: X
     &out
      real(8), intent(in), dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: Y
     &out
      real(8), intent(out), dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: 
     &Fout
        end subroutine Sub_Loop_cinterp2d

      end interface
      integer(4), intent(in) :: LBx
      integer(4), intent(in) :: UBx
      integer(4), intent(in) :: LBy
      integer(4), intent(in) :: UBy
      integer(4), intent(in) :: Istr
      integer(4), intent(in) :: Iend
      integer(4), intent(in) :: Jstr
      integer(4), intent(in) :: Jend
      real(8), intent(in), dimension(LBx:UBx) :: Xinp
      real(8), intent(in), dimension(LBy:UBy) :: Yinp
      real(8), intent(in), dimension(LBx:UBx,LBy:UBy) :: Finp
      real(8), intent(in), dimension(-1: Agrif_tabvars_i(2)%iarray0+2+ A
     &grif_tabvars_i(13)%iarray0,-1: Agrif_tabvars_i(4)%iarray0+2+ Agrif
     &_tabvars_i(14)%iarray0) :: Xout
      real(8), intent(in), dimension(-1: Agrif_tabvars_i(2)%iarray0+2+ A
     &grif_tabvars_i(13)%iarray0,-1: Agrif_tabvars_i(4)%iarray0+2+ Agrif
     &_tabvars_i(14)%iarray0) :: Yout
      real(8), intent(out), dimension(-1: Agrif_tabvars_i(2)%iarray0+2+ 
     &Agrif_tabvars_i(13)%iarray0,-1: Agrif_tabvars_i(4)%iarray0+2+ Agri
     &f_tabvars_i(14)%iarray0) :: Fout
      

        call Sub_Loop_cinterp2d(LBx,UBx,LBy,UBy,Xinp,Yinp,Finp,Istr,Iend
     &,Jstr,Jend,Xout,Yout,Fout, Agrif_tabvars_i(14)%iarray0, Agrif_tabv
     &ars_i(4)%iarray0, Agrif_tabvars_i(13)%iarray0, Agrif_tabvars_i(2)%
     &iarray0)

      end


      subroutine Sub_Loop_cinterp2d(LBx,UBx,LBy,UBy,Xinp,Yinp,Finp,Istr,
     &Iend,Jstr,Jend,Xout,Yout,Fout,padd_E,Mm,padd_X,Lm)

      use Agrif_Types, only : Agrif_tabvars

      implicit none

      integer(4), parameter :: LLm0 = 41
      integer(4), parameter :: MMm0 = 42
      integer(4), parameter :: N = 32
      integer(4), parameter :: NP_XI = 1
      integer(4), parameter :: NP_ETA = 4
      integer(4), parameter :: NNODES = NP_XI*NP_ETA
      integer(4), parameter :: NPP = 1
      integer(4), parameter :: NSUB_X = 1
      integer(4), parameter :: NSUB_E = 1
      integer(4), parameter :: NWEIGHT = 1000
      integer(4), parameter :: stdout = 6
      integer(4), parameter :: Np = N+1
      integer(4), parameter :: NSA = 28
      real, parameter :: Vtransform = 2
      integer(4), parameter :: itemp = 1
      integer(4), parameter :: ntrc_temp = 1
      integer(4), parameter :: ntrc_salt = 1
      integer(4), parameter :: ntrc_pas = 0
      integer(4), parameter :: ntrc_bio = 0
      integer(4), parameter :: ntrc_subs = 0
      integer(4), parameter :: ntrc_substot = 0
      integer(4), parameter :: ntrc_sed = 0
      integer(4), parameter :: NTA = itemp+ntrc_salt
      integer(4), parameter :: NT = itemp+ntrc_salt+ntrc_pas+ntrc_bio+nt
     &rc_sed
      integer(4), parameter :: NTot = NT
      integer(4), parameter :: isalt = itemp+1
      integer(4), parameter :: ntrc_diabio = 0
      integer(4), parameter :: ntrc_diats = 0
      integer(4), parameter :: ntrc_diauv = 0
      integer(4), parameter :: ntrc_diavrt = 0
      integer(4), parameter :: ntrc_diaek = 0
      integer(4), parameter :: ntrc_diapv = 0
      integer(4), parameter :: ntrc_diaeddy = 0
      integer(4), parameter :: ntrc_surf = 0
      real(8), parameter :: C01 = 1.0D0/48.0D0
      real(8), parameter :: C02 = 1.0D0/32.0D0
      real(8), parameter :: C03 = 0.0625D0
      real(8), parameter :: C04 = 1.0D0/6.0D0
      real(8), parameter :: C05 = 0.25D0
      real(8), parameter :: C06 = 0.5D0
      real(8), parameter :: C07 = 0.3125D0
      real(8), parameter :: C08 = 0.625D0
      real(8), parameter :: C09 = 1.5D0
      real(8), parameter :: C10 = 13.0D0/24.0D0
      real(8), parameter :: LIMTR = 3.0D0
      real(8), parameter :: spv = 0.0D0
      integer(4) :: padd_E
      integer(4) :: Mm
      integer(4) :: padd_X
      integer(4) :: Lm
      integer(4), intent(in) :: LBx
      integer(4), intent(in) :: UBx
      integer(4), intent(in) :: LBy
      integer(4), intent(in) :: UBy
      integer(4), intent(in) :: Istr
      integer(4), intent(in) :: Iend
      integer(4), intent(in) :: Jstr
      integer(4), intent(in) :: Jend
      real(8), intent(in), dimension(LBx:UBx) :: Xinp
      real(8), intent(in), dimension(LBy:UBy) :: Yinp
      real(8), intent(in), dimension(LBx:UBx,LBy:UBy) :: Finp
      real(8), intent(in), dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: X
     &out
      real(8), intent(in), dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: Y
     &out
      real(8), intent(out), dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: 
     &Fout

                   

                                            
                             
                                            

                                                             
                                                            
                                           
                                                               
                                    
                                      
                                                          

                       

                                    

                        
                              

                                          
                                             
                          

                        

                                          
                        

                                                        
                      
                              

                                       
                                                                     
                                          
                         

                             

                             

                            

                            

                                             

                            

                                     

                                                               

                         

                                                      
                                                      
                                                                  
                               

                               

                              

                              

                               

                              

                              

                                

                             

                                            
                                               
                                            
                                               
                                                                 
                                                 
                                 
                                               
                                               
                                               
                                               
                                               
                                               
                                             
                                             
                                             
                                             
                                             
                                               
                                               
                                               
                                               
                                               
                                               
                                               
                                               
                                               
                                               
                                               
                                               
                                              
                                                  
                                                  
                                                  
                                                  
                                                 
                                                 
                                               
                                               
                                                      
                                                 
                                                 
                                                 
                                                 
                                                 
                                                 
                                                 
                                                           
                                                           
                                                           
                                    
                                                
                                                
                                                
                                                
                                                 
                           
                           
                           
                           
                             
                                              
                          
                                                  
                                                      
                                                
                                                
                                                        
                                                                      
                                                                      
                                                                       
                                                           
      integer(4) :: i
      integer(4) :: ic
      integer(4) :: iter
      integer(4) :: i1
      integer(4) :: i2
      integer(4) :: j
      integer(4) :: jc
      integer(4) :: j1
      integer(4) :: j2
      integer(4) :: ii
      integer(4) :: jj
                                         
      real(8) :: a11
      real(8) :: a12
      real(8) :: a21
      real(8) :: a22
                                         
      real(8) :: e11
      real(8) :: e12
      real(8) :: e21
      real(8) :: e22
                                                                
      real(8) :: cff
      real(8) :: d1
      real(8) :: d2
      real(8) :: dfc
      real(8) :: dx
      real(8) :: dy
      real(8) :: eta
      real(8) :: xi
      real(8) :: xy
      real(8) :: yx
                                                                       
                                                                 
 
      real(8) :: f0
      real(8) :: fx
      real(8) :: fxx
      real(8) :: fxxx
      real(8) :: fxxy
      real(8) :: fxy
      real(8) :: fxyy
      real(8) :: fy
      real(8) :: fyy
      real(8) :: fyyy
                                                    
                                                    
                                                
                                                   
                                              
                                             
                                                
                                               
                                             
                                                     
                                               
                                             
                                                         
      real(8), dimension(-1:2,-1:2) :: dfx
      real(8), dimension(-1:2,-1:2) :: dfy
      real(8), dimension(-1:2,-1:2) :: ff
      DO j=Jstr,Jend
        DO i=Istr,Iend
           DO ii=LBx,(UBx-1)
             if ((Xinp(ii).le.Xout(i,j)).and.
     &           (Xinp(ii+1).gt.Xout(i,j))) then
               i1=ii
               i2=ii+1
               goto 10
             endif
           enddo
           print*, 'Did not find i1 and i2',
     &           Istr,Iend,Jstr,Jend,i,j,Xout(i,j),Xout(i-1,j)
           goto 100
 10        continue
           DO jj=LBy,UBy-1
             if ((Yinp(jj).le.Yout(i,j)).and.
     &           (Yinp(jj+1).gt.Yout(i,j))) then
               j1=jj
               j2=jj+1
               goto 20
             endif
           enddo
           print*, 'Did not find j1 and j2'
           goto 100
 20        continue
          IF (((LBx.le.i1).and.(i1.le.UBx)).and.
     &        ((LBy.le.j1).and.(j1.le.UBy))) THEN
            xy=Xinp(i2)-Xinp(i1)-Xinp(i2)+Xinp(i1)
            yx=Yinp(j2)-Yinp(j2)-Yinp(j1)+Yinp(j1)
            dx=Xout(i,j)-0.25D0*(Xinp(i2)+Xinp(i1)+
     &                         Xinp(i2)+Xinp(i1))
            dy=Yout(i,j)-0.25D0*(Yinp(j2)+Yinp(j2)+
     &                         Yinp(j1)+Yinp(j1))
            e11=0.5D0*(Xinp(i2)-Xinp(i1)+Xinp(i2)-Xinp(i1))
            e12=0.5D0*(Xinp(i2)+Xinp(i1)-Xinp(i2)-Xinp(i1))
            e21=0.5D0*(Yinp(j2)-Yinp(j2)+Yinp(j1)-Yinp(j1))
            e22=0.5D0*(Yinp(j2)+Yinp(j2)-Yinp(j1)-Yinp(j1))
            cff=1.0D0/(e11*e22-e12*e21)
            xi=cff*(e22*dx-e12*dy)
            eta=cff*(e11*dy-e21*dx)
            DO iter=1,4
              d1=dx-e11*xi-e12*eta-xy*xi*eta
              d2=dy-e21*xi-e22*eta-yx*xi*eta
              a11=e11+xy*eta
              a12=e12+xy*xi
              a21=e21+yx*eta
              a22=e22+yx*xi
              cff=1.0D0/(a11*a22-a12*a21)
              xi =xi +cff*(a22*d1-a12*d2)
              eta=eta+cff*(a11*d2-a21*d1)
            END DO
            DO jc=-1,2
              DO ic=-1,2
                ff(ic,jc)=Finp(MAX(1,MIN(UBx,i1+ic)),
     &                         MAX(1,MIN(UBy,j1+jc)))
              END DO
            END DO
            f0=C07*(ff(1,1)+ff(1,0)+ff(0,1)+ff(0,0))-
     &         C02*(ff(2,0)+ff(2,1)+ff(1,2)+ff(0,2)+
     &              ff(-1,1)+ff(-1,0)+ff(0,-1)+ff(1,-1))
            fx=C08*(ff(1,1)+ff(1,0)-ff(0,1)-ff(0,0))-
     &         C01*(ff(2,1)+ff(2,0)-ff(-1,1)-ff(-1,0))-
     &         C03*(ff(1,2)-ff(0,2)+ff(1,-1)-ff(0,-1))
            fy=C08*(ff(1,1)-ff(1,0)+ff(0,1)-ff(0,0))-
     &         C01*(ff(1,2)+ff(0,2)-ff(1,-1)-ff(0,-1))-
     &         C03*(ff(2,1)-ff(2,0)+ff(-1,1)-ff(-1,0))
            fxy=ff(1,1)-ff(1,0)-ff(0,1)+ff(0,0)
            fxx=C05*(ff(2,1)-ff(1,1)-ff(0,1)+ff(-1,1)+
     &               ff(2,0)-ff(1,0)-ff(0,0)+ff(-1,0))
            fyy=C05*(ff(1,2)-ff(1,1)-ff(1,0)+ff(1,-1)+
     &               ff(0,2)-ff(0,1)-ff(0,0)+ff(0,-1))
            fxxx=C06*(ff(2,1)+ff(2,0)-ff(-1,1)-ff(-1,0))-
     &           C09*(ff(1,1)+ff(1,0)-ff(0,1)-ff(0,0))
            fyyy=C06*(ff(1,2)+ff(0,2)-ff(1,-1)-ff(0,-1))-
     &           C09*(ff(1,1)-ff(1,0)+ff(0,1)-ff(0,0))
            fxxy=C06*(ff(2,1)-ff(1,1)-ff(0,1)+ff(-1,1)-
     &                ff(2,0)+ff(1,0)+ff(0,0)-ff(-1,0))
            fxyy=C06*(ff(1,2)-ff(1,1)-ff(1,0)+ff(1,-1)-
     &                ff(0,2)+ff(0,1)+ff(0,0)-ff(0,-1))
            Fout(i,j)=f0+
     &                fx*xi+
     &                fy*eta+
     &                C06*fxx*xi*xi+
     &                fxy*xi*eta+
     &                C06*fyy*eta*eta+
     &                C04*fxxx*xi*xi*xi+
     &                C06*fxxy*xi*xi*eta+
     &                C04*fyyy*eta*eta*eta+
     &                C06*fxyy*xi*eta*eta
          END IF
        END DO
      END DO
      RETURN
 100  continue
      print*, 'error in cinterp2d'
             

      end subroutine Sub_Loop_cinterp2d


