CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCC
CCCCCCCCCCC  PROGRAMA PARA EXTRAIR VARIAS SECOES DE VELOCIDADES
CCCCCCCCCCC  SIMULTANEAMENTE
CCCCCCCCCCC  (BINARIO - ACESSO DIRETO), E  SALVAR EM UM BINARIO ([.a])
CCCCCCCCCCC
CCCCCCCCCCC  MARIELA GABIOUX (SETEMBRO DE 2011)
CCCCCCCCCCC  Baseado na leitura/escrita do hycom (mod_za.f)
CCCCCCCCCCC  e nos scripts criados pelo Prof. Afonso Paiva e Vladimir
CCCCCCCCCCC  Santos da Silva
CCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCC
CCCCCCCCCCC  Versao: tira_secao_interp_vsc_mg.f
CCCCCCCCCCC     le arquivo para varios tempos, extrai uma secao de veloc
CCCCCCCCCCC     zonal interpolada na vertical e salva em 1 arquivo de saida
CCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      program calcula_secao_zonal_VTOT_archv_ATLa008

      USE IFPORT
      implicit none

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c ... DEFINE INDICE DO REGISTRO A SER LIDO NO ARQUIVO DO HYCOM

c ... MODELO 1/12 novo atlantico (versao 2.2.18 do HYCOM)
      integer(4), parameter:: JDM=1100, IDM=1500, KDM=41
      integer(4):: MONTG_INDEX=0, SSH_INDEX=1
      integer(4):: ONETA_INDEX=2, SURFLX_INDEX=3
      integer(4):: WTRFLX_INDEX=4, SALFLX_INDEX=5
      integer(4):: BLDEPTH_INDEX=6, MLDEPTH_INDEX=7

c      integer(4):: pk14=6

      integer(4):: UBARO_INDEX=11, VBARO_INDEX=12
      integer(4):: UVEL_INDEX=13, VVEL_INDEX=14, LT_INDEX=15
      integer(4):: TEMP_INDEX=16, SALT_INDEX=17

      integer(4):: NUM3DF=5
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      integer(4), parameter:: n2drec=((idm*jdm+4095)/4096)*4096
      integer(4):: k, i, j, INDEX, ios, recnum, ind
c ----------------------------------------------------------
c ---- define parametros da secao a extrair CB
      integer(4), parameter:: KK=33 !185 !276  !33
      integer(4), parameter:: i1=1  !-36W 
      integer(4), parameter:: i2=IDM  !-20W

      integer(4), parameter:: iii=i2-i1+1
      integer(4), parameter:: n2drec2=((iii*kk+4095)/4096)*4096

c -----------------------------------------------------------
c      integer(4), parameter:: nrolatmax=11
      integer(4), parameter:: nrolatmax=1 !5S e 10S
      integer(4):: nn,jsec,jsec0(nrolatmax),n1,n2
      character expt0*10,expt1*23

      integer(4):: archv_idin=12, archv_idout=13, archv_log=11
      integer(4)::archv_b=16
      integer(4):: IJDM=IDM*JDM 
      integer(4):: ano,dia,anoini,anoend,diaini,diaend,inter
      character arq_in*150, arq_out*150, arq_log*150,work*15
      character arq_b*150
      real(4):: h(idm,jdm),lt(idm,jdm),h1(iii,kdm),lt1(iii,kdm)
      real(4):: h0(idm,jdm),lat(nrolatmax)
      real(4):: h2(kdm),lt2(kdm),h_interp1(kk),h_interp2(iii,kk)
      real(4):: prof(kk)
      real(4):: spval=1.e20
      real(4):: bb=1/2.0001

c -----------------------------------------------------------
c === niveis levitus
      data (prof(k),k=1,kk) 
     &/0,10,20,30,50,75,100,125,150,200,250,300,400,500,600,700,
     &800,900,1000,1100,1200,1300,1400,1500,1750,2000,2500,3000,3500,
     &4000,4500,5000,5500/ !profundidades Levitus
 
c =============================================
      !expt0='ATL12'

      expt1='GLB025t09'

        anoini=2006
        anoend=2007

C        diaini=1
        inter=1

c --- ######################################
c       lat(1)=5   ! 5S
       lat(1)=-67   ! 10S
c --- ######################################
c       jsec0(1)=757  ! indice de latitude 5 S
       jsec0(1)=112 !indice de latitude S 10

c --- #####################################
c --- define seção inicial e final
      n1=1
      n2=nrolatmax
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!! inicia loop latitudes	 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c      do 15 nn = 1,nrolatmax
      do 15 nn = n1,n2
        jsec=jsec0(nn)

      do 20 ano = anoini,anoend

        IF (ano .EQ. 2004 .OR. ano .EQ. 2016) THEN
           diaini = 1
           diaend=366   
	ELSE
           diaini = 1
           diaend=365
        ENDIF

        recnum=0

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c ... abre log de saida

       write(arq_log,'("seczon_j_",i4.4,"_",a9,"_V_",i4.4,"_",
     &i3.3,"a",i3.3,".log")')jsec,expt1,ano,diaini,diaend

      open(archv_log,file=arq_log,status='unknown')


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c ... abre arquivo de saida ascii .b

       write(arq_b,'("seczon_j_",i4.4,"_",a9,"_V_",i4.4,"_",
     &i3.3,"a",i3.3,".b")')jsec,expt1,ano,diaini,diaend

      open(archv_b,file=arq_b,status='unknown')

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      write(archv_b,*)'*************************** '
      write(archv_b,*)' '
      write(archv_b,*)' Informacoes do arquivo'
      write(archv_b,*)' '
      write(archv_b,102)lat(nn),jsec
      write(archv_b,*)' '
      write(archv_b,103)ano,diaini,diaend
      write(archv_b,*)' '
      write(archv_b,104)inter
      write(archv_b,*)' '
      write(archv_b,*)'Interpolada na grade de p'
      write(archv_b,*)' '
      write(archv_b,105)kk
      write(archv_b,*)' '
      write(archv_b,*)'*************************** '
102   format('  Seção zonal de velocidade latitude: ',f7.2,
     &' indice ',i4.4)
103   format('  Ano: ',i4,' dia inicial: ',i3.3,' dia final: ',i3.3)
104   format('  Informação salva a cada ',i2.2,' dias')
105   format('  total de profundidades variaveis: ',i3.3)


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c ... abre arquivo para salvar resultados (binario)

      write(arq_out,'("seczon_j_",i4.4,"_",a9,"_V_",i4.4,"_",i3.3,
     &"a",i3.3,".a")')jsec,expt1,ano,diaini,diaend


      open(archv_idout,file=arq_out,status='unknown',form='unformatted'
     &,recl=n2drec2,access='DIRECT',action='WRITE')
      write(*,'("resultados em: ",a)')arq_out

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!! inicia loop temporal       (ex: todos os dias de 2006)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      do 30 dia = diaini,diaend,inter
        recnum=recnum+1 ! No do registro a ser salvo

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!! Experimento 1/16 sem assimilação 2005 !!!!!!!!!!!!!!!!!!!!!!!
! Base Hidrodinâmica REMO 3
      write(arq_in,'("
     &/pesq/share/monan/ocean_model/hycom/GLBa025/dados/
     &data_teste09_",i4.4,"/
     &archv.",i4.4,"_",i3.3,"_00.a")')ano,ano,dia


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c ===========================================================
c ===========================================================

c ... abre arquivo
      open(unit=archv_idin,file=arq_in,status='old',form='unformatted',
     &action='READ',recl=n2drec,access='direct',iostat=ios) 

c      write(archv_log,*)'ano,dia = ',ano,dia
      write(archv_log,'(a)')arq_in
c      write(*,'(a)')arq_in

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c ... leitura dos campos 2D

      INDEX = VBARO_INDEX
      call read2d(h0,archv_idin,idm,jdm,n2drec,index)

C ... define o campo a ser lido
      do 10 k=1,kdm
      INDEX = ((k-1)*NUM3DF+VVEL_INDEX)

c ... le registro - vetor a(n2drec)
      call read2d(h,archv_idin,idm,jdm,n2drec,index)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c ... exemplo: modifica spval de terra, mantem TSM como esta
c     mas se for SSH ou ML_depth, ou dp, entao converte unidades
c ... O importante é que neste ponto do programa voce tem uma matriz 
c     h(idm,jdm) contendo o campo 2D lido. Faca o que quiser com ela.

      do i=i1,i2

        if(h0(i,jsec).gt.1.e10) then
          h1(i-i1+1,k) = spval
        else
          h1(i-i1+1,k) = ( h0(i,jsec) + h(i,jsec) + h0(i,jsec+1) + 
     & h(i,jsec+1))*0.5       
        endif
      enddo

C ... define o campo a ser lido
      INDEX = ((k-1)*NUM3DF+LT_INDEX)
  
c ... le registro - vetor a(n2drec)
      call read2d(lt,archv_idin,idm,jdm,n2drec,index)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      do i=i1,i2
        if(lt(i,jsec).gt.1.e10) then
          lt1(i-i1+1,k) = spval
        else
          lt1(i-i1+1,k) = lt(i,jsec)  / 9806.          
        endif
      enddo

 10   continue

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c ... fazendo a interpolacao vertical para cada ponto
c ... entrada: nome do arquivo de saida 
c ... campo 1D (pontual) de espessura de camada (lt2);
c ... variavel do ponto em todas as profundidades (h2);
c ... profundidades a ser interpolado (prof);
c ... numero total de camadas (KDM)
c ... numero de profundidades a ser interpoladas
c ... bb= coef para dividir a camada, bb deve ser menor a 1/2 pois nesse caso
c ... os dois pontos estariam no centro da camada.
c ... Sugestao bb=1/(2.0001)
c       write(*,*)'interpola'
      do i=1,iii
      do k=1,kdm
        h2(k) = h1(i,k)
        lt2(k) = lt1(i,k)
      enddo

      call interpz(h_interp1,lt2,h2,prof,kdm,KK,bb,i)
      
      do k=1,kk
      h_interp2(i,k) = h_interp1(k);
      enddo
      enddo ! do de j

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c     salva arquivo binario com campo 2D extraido/manipulado
      call  write2d(h_interp2,archv_idout,iii,kk,n2drec2,recnum)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c ... fecha arquivo de leitura antes de ir para proximo tempo
      close(archv_idin)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!! fecha loop de tempo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
30    continue
c ... fecha arquivo de saida, log e b
      close(archv_idout)
      close(archv_log)
      close(archv_b)
20    continue
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
15    continue
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c ... encerra programa
      stop
      end   
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!   subrotina read2d
!!!!!!!!!!   le campo 2D (camada) em vetor a(n2drec), convert endian
!!!!!!!!!!   e transforma vetor a(n2drec) em matriz aa(idm,jdm),
!!!!!!!!!!   retornando [aa] para o programa principal
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      subroutine read2d(aa,archv_idin,idm,jdm,n2drec,index)

      implicit none
      integer(4):: i,j,idm,jdm,n2drec,index,ios
      integer(4):: archv_idin
      real(4):: a(n2drec),aa(idm,jdm)

c ... le registro - vetor a(n2drec)
      read(archv_idin,rec=INDEX+1,iostat=ios)a
      call zaio_endian(a,n2drec)
c ... transforma vetor lido em matriz 2D
      do j=1,jdm
      do i=1,idm
        aa(i,j) = a(i+(j-1)*idm)
      enddo
      enddo
      return
      end subroutine read2d

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!   subrotina write2d
!!!!!!!!!!   recebe campo 2D (camada) na matriz aa(idm,jdm), corrige
!!!!!!!!!!   endian, transforma no vetor a(n2drec) e salva este vetor
!!!!!!!!!!   no registro recnum, ou seja, salva campo 2D
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      subroutine write2d(aa,archv_idout,idm,jdm,n2drec,recnum)

      implicit none
      integer(4):: i,j,idm,jdm,n2drec,recnum,ios
      integer(4):: archv_idout
      real(4):: a(n2drec),aa(idm,jdm)

c ... transforma campo 2D em vetor a(n2drec) para salvar
      do j=1,jdm
      do i=1,idm
        a(i+(j-1)*idm) = aa(i,j)
      enddo
      enddo
c ... corrige endian e salva arquivo binario
      call zaio_endian(a,n2drec)
      write(unit=archv_idout,rec=recnum,iostat=ios)a
      return
      end subroutine write2d

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!
!!!!!!!!!   zaio_endian(a,n)
!!!!!!!!!   - swap the endian-ness of the array.
!!!!!!!!!   - assumes integer(kind=1) and integer(kind=4) ocupy
!!!!!!!!!     one and four bytes respectively.
!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      subroutine zaio_endian(a,n)
      implicit none
      integer,         intent(in)    :: n
      integer(kind=4), intent(inout) :: a(n)  ! 4-bytes
      integer         k
      integer(kind=4) ii4,   io4     ! 4-bytes
      integer(kind=1) ii1(4),io1(4)  ! 1-byte
      equivalence    (ii4,ii1(1)), (io4,io1(1))  ! non-standard f90
      do k= 1,n
        ii4 = a(k)
        io1(1) = ii1(4)
        io1(2) = ii1(3)
        io1(3) = ii1(2)
        io1(4) = ii1(1)
        a(k) = io4
      enddo
      return
      end subroutine zaio_endian
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	subroutine interpz(dadoi,hm,dado,zz,Kt,M,bb,i)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c...function [dadoi]=interpz(hm,dado,zz,K,M,bb);
c...
c... Interpolacao dos dados de camadas para niveis
c...
c... Variaveis necessarias:
c... hm = espessuras das camadas
c... dado = variavel a ser interpolada
c... M = n. de profundidades na saida
c... Kt = total de camadas
c... zz = novo eixo de profundidades;
c... bb= coef para dividir a camada, bb deve ser menor a 1/2 pois nesse caso
c... os dois pontos estariam no centro da camada.
c... Sugestao bb=1/(2.0001)


      implicit none
      integer(4):: mm,M,kk,kt,i,j,n,k
      real(4):: bb
      real(4):: hh(kt),hm(kt),p(kt+1),z(2*kt+1),dd(2*kt+1)
      real(4):: dado(kt),zz(M)
      real(4):: dadoi(M)

c ... retira valores extremos

      do kk=2,kt
         if (dado(kk).gt..5e2 .or. dado(kk).lt.-.5e2)
     &   dado(kk)=dado(kk-1)  
      enddo
c...Calcula profundidades acumuladas 
      p(1)=0.0
    
      do kk=1,Kt
	 hh(kk)=max( hm(kk),0.1 )
	 p(kk+1)=p(kk)+hh(kk)  
      enddo

c...localizando prof dos resultados e interpolando
       
	 n=1 
         z(n)=0.0 
         dd(1)=dado(1)     

	 do kk=1,Kt
	    n=n+1;
	    z(n)=p(kk)+hh(kk)*bb
	    dd(n)=dado(kk)
	    n=n+1
	    z(n)=p(kk+1)-hh(kk)*bb
	    dd(n)=dado(kk)
	 enddo               
                    
	 z(n+1)=p(kt+1) 
         dd(n+1)=dd(n)
            
	 if (zz(M) .gt. z(2*kt+1)) then
         z(2*kt+1)=zz(M)+1
         endif

c	do k=1,2*kt+1
c         if (i.eq.361) write (101,*)k,dd(k),z(k)
c        enddo


	 call interp_1D(dadoi,z,dd,zz,kt,M)  

c        do k=1,33
c         if (i.eq.361) write (102,*)zz(k),dadoi(k)
c        enddo
	
      return
      end subroutine interpz


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine interp_1D(aa,z,dd,zz,kt,M)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc


      implicit none
      integer(4):: i,j,kt,M,diff
      real(4):: bb,cc
      real(4):: aa(M),dd(2*kt+1),zz(M),z(2*kt+1)

      
      do i=1,M
         j=1
         diff = zz(i) - z(j)
         do while (diff .ge. 0 .and. j .le. 2*kt)
         j=j+1
         diff = zz(i) - z(j)
         enddo

         if (j.eq.1) then
         aa(i) = dd(j)
         elseif (j.eq. 2*kt+1 .and. diff.ge.0) then
         aa(i) = dd(j)   
         else   
         bb = abs(zz(i) - z(j-1))
         cc = abs(zz(i) - z(j))
         aa(i) = (dd(j-1)*cc + dd(j)*bb)/(bb+cc)    
         endif     
      enddo

      return
      end subroutine interp_1D
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
