CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCC
CCCCCCCCCCC  PROGRAMA PARA LER RESULTADOS DO HYCOM (BINARIO - ACESSO 
CCCCCCCCCCC  DIRETO), CALCULAR ALGUMA COISA OU EXTRAIR ALGUM CAMPO E
CCCCCCCCCCC  SALVAR EM UM BINARIO (O FAMOSO [.a]) 
CCCCCCCCCCC
CCCCCCCCCCC  AFONSO PAIVA (SETEMBRO DE 2010)
CCCCCCCCCCC  Baseado na leitura/escrita do hycom (mod_za.f) 
CCCCCCCCCCC  Muchas Gracias a Guedes e Vladimir pela ajuda
CCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCC
CCCCCCCCCCC  Este é somente um programa básico para ser adaptado a
CCCCCCCCCCC  cada necessidade de processamento dos resultados
CCCCCCCCCCC
CCCCCCCCCCC  Funciona tanto no cluster como aqui. A idéia basica é:
CCCCCCCCCCC    - facilidade de extrair campo para analise (posterior no 
CCCCCCCCCCC      matlab ou python) sem ter que ficar fazendo mil copias 
CCCCCCCCCCC      dos arquivos de resultados em todos os computadores  
CCCCCCCCCCC    - provavelmente o fortran é mais eficiente do que o
CCCCCCCCCCC      matlab no tratamento de grandes arquivos, enquanto o
CCCCCCCCCCC      matlab será mais util para analise e plotagem
CCCCCCCCCCC    - facilidade de calcular campos medios, transporte em 
CCCCCCCCCCC      secoes, interpolacao para niveis, etc
CCCCCCCCCCC    - o mesmo procedimento é facilmente adaptavel para gerar
CCCCCCCCCCC      forcantes a partir de NCEP ou outras bases, sem ficar 
CCCCCCCCCCC      de todos aqueles programas complicados do hycom
CCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCC
CCCCCCCCCCC  Versao: tira_prof_especifica_U
CCCCCCCCCCC  Como va1, mas:
CCCCCCCCCCC     le arquivo para varios tempos, extrai um campo 2d
CCCCCCCCCCC     de cada um (ex: TSM), cria matriz h(idm,jdm) e calcula
CCCCCCCCCCC     algo com ela, e salva campos 2D em 1 arquivo de saida
CCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      program tira_prof_especifica_U

      USE IFPORT
      implicit none

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c ... DEFINE INDICE DO REGISTRO A SER LIDO NO ARQUIVO DO HYCOM
      integer(4), parameter:: JDM=1100, IDM=1500, KDM=32
      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
      integer(4):: COVICE_INDEX=8, THKICE_INDEX=9
      integer(4):: TEMICE=10
      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 !NUmber of 3D fields set = 0 if you want
                             ! to read only the 2d fields

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      integer(4), parameter:: n2drec=((idm*jdm+4095)/4096)*4096
      integer(4):: k, i, j, INDEX, ios, recnum, ind, index0
      integer(4), parameter:: kk=1
      integer(4):: archv_idin=12,archv_idout=13,archv_log=11,archv_b=16
      integer(4):: IJDM=IDM*JDM 
      integer(4):: ano,dia,anoini,anoend,diaini,diaend,prof_name
      integer(4):: inter,mes,dd,diad,dini,dend,hora
      integer(4):: bi !! 1 para ano bissexto e 0 caso contrario
      character arq_in*200, arq_out*150,arq_b*120
      character arq_log*150
      character expt1*23
      real(4):: h(idm,jdm),lt(idm,jdm),h0(idm,jdm)
      real(4):: h1(kdm,idm,jdm),lt1(kdm,idm,jdm)
      real(4):: h2(kdm),lt2(kdm),h_interp1(kk)
      real(4):: h_interp2(idm,jdm)

      real(4):: profs=0  !prof onde será interpolada
c      real(4):: profs=15  !prof onde será interpolada
c      real(4):: profs=50
c      real(4):: profs=300 
c      real(4):: profs=1500  !prof onde será interpolada
c      real(4):: profs=2000  !prof onde será interpolada
c      real(4):: profs=4200  !prof onde será interpolada

      real(4):: spval=1.e20
      real(4):: bb=1/2.0001

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c --- define parametros

      expt1='GLb25'
      prof_name=int(profs)
      inter=1
      anoini=2005
      anoend=2005
      
      diaini=1
c      diaend=365
      inter=1
c      recnum=0
      do ano = anoini,anoend,1
      recnum=0  
c      do dia = diaini,diaend,1
c        recnum=recnum+1
        IF (ano .EQ. 2004 .OR. ano .EQ. 2008 .OR. ano .EQ. 2012 .OR.ano
     &.EQ. 2016 ) THEN
           diaend=366
        ELSEIF (ano .EQ. 2025) THEN
           diaend=238
        ELSE
           diaend=365
c           diaend=30 !3
       ENDIF
!      do dia = diaini,diaend,1
!       recnum=recnum+1

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      write(arq_log,'("Vinterp_",a5,"_",i4.4,"m_",
     &i4.4,".log")')expt1,prof_name,ano
      open(archv_log,file=arq_log,status='unknown')

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

       write(arq_b,'("Vinterp_",a5,"_",i4.4,"m_",
     &i4.4,".b")')expt1,prof_name,ano
      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,*)''
      write(archv_b,*)'*************************** '
      write(archv_b,102)profs
      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,*)'*************************** '

102   format('  Componente V interpolada na Profundidade [m]: ',f7.2)
103   format('  Data:',i4,' dia inicial: ',i3.3,' dia final: ',i3.3)
104   format('  Campos salvos a cada ',i2.2,' dias')

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

       write(arq_out,'("Vinterp_",a5,"_",i4.4,"m_",
     &i4.4,".a")')expt1,prof_name,ano
      open(archv_idout,file=arq_out,status='unknown',form='unformatted'
     &,recl=n2drec,access='DIRECT',action='WRITE')
      write(*,'("resultados em: ",a)')arq_out

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!! inicia loop temporal       (ex: todos os dias de 2006)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c        recnum=0
c      do dia = diaini,diaend,inter
c        recnum=recnum+1 ! No do registro a ser salvo
c        write(*,*)recnum
      do dia = diaini,diaend,1
c       do hora=0,12,12
       recnum=recnum+1
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c ... monta nome e abre arquivo de dados
      write(arq_in,'("/mnt/beegfs/efraime.daniel/hycom/GLBa025
     &/teste19_ERA5_Newbatim/data_",i4.4,"/
     &archv.",i4.4,"_",i3.3,"_00.a")')ano,dia
c ... abre arquivo
      open(unit=archv_idin,file=arq_in,status='old',form='unformatted',
     &action='READ',recl=n2drec,access='direct',iostat=ios)

      write(archv_log,'(a)')arq_in

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

c ... le registro - vetor a(n2drec)
      INDEX0 = VBARO_INDEX
c      INDEX0 = UBARO_INDEX
      call read2d(h0,archv_idin,idm,jdm,n2drec,index0)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
C ... define o campo a ser lido
      do k=1,kdm
      INDEX = ((k-1)*NUM3DF+VVEL_INDEX)
c      INDEX = ((k-1)*NUM3DF+UVEL_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=1,idm-1
      do j=1,jdm-1

        if(h(i,j).gt.1.e10) then
          h1(k,i,j) = spval
        else
           h1(k,i,j) = ( h(i,j) + h0(i,j) +  h(i,j+1) + h0(i,j+1) )*0.5  ! p/ interpolar V       
c          h1(k,i,j) = ( h(i,j) + h0(i,j) +  h(i+1,j) + h0(i+1,j) )*0.5  ! p/ interpolar U
        endif

      enddo
      enddo
      
      do j=1,jdm
      h1(k,idm,j) = spval
      enddo

      do i=1,idm
      h1(k,i,jdm) = spval
      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=1,idm
      do j=1,jdm

        if(lt(i,j).gt.1.e10) then
          lt1(k,i,j) = spval
        else
          lt1(k,i,j) = lt(i,j)  / 9806.          
        endif

      enddo
      enddo

c 10   continue
      enddo !camada
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
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)
    
      do i=1,idm
      do j=1,jdm
      do k=1,kdm

        h2(k) = h1(k,i,j)
        lt2(k) = lt1(k,i,j)

      enddo

      call interpz(h_interp1,lt2,h2,profs,kdm,1,bb,i)
      
      do k=1,1,1
c      h_interp2(k,i,j) = h_interp1(k);
       h_interp2(i,j) = h_interp1(k);
      enddo

      enddo
      enddo

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c     salva arquivo binario com campo 2D extraido/manipulado      
        write(archv_log,*)'registro para escrita (recnum) = ',recnum
        call  write2d(h_interp2,archv_idout,idm,jdm,n2drec,recnum)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c ... fecha arquivo de leitura antes de ir para proximo tempo
      close(archv_idin)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!! fecha loop de tempo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!c30    continue
c      enddo !hora
      enddo !dia
c ... fecha arquivo de saida e log
c      close(archv_idout)
c      close(archv_log)
c      close(archv_b)

c      enddo    
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c ... fecha arquivo de saida e log
      close(archv_idout)
      close(archv_log)
      close(archv_b)
      enddo !ano
 
!c        20    c
!c ontinue

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
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
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!