!*******************************************************************************
!*                                 MGRADS									   *
!*																			   *
!* Modulo de funcoes e subrotinas uteis para gravacao e leitura de dados       *
!*                    meteorologicos no formato do GRADS                       *
!*                                                                             *                                                                             *
!*                                                                        	   *
!*      Copyright (C) 2005 	Sergio Henrique S. Ferreira                        *
!* 																			   *
!*      MCT-INPE-CPTEC-Cachoeira Paulista, Brasil                              *
!*                                                                             *
!*-----------------------------------------------------------------------------*
!*                                                                             *
!*   SUB-ROTINAS PUBLICAS                                                      *
!*                                                                             *
!*     SAVEBIN - grava linhas de dados observacionais em formato binario do    * 
!*               grads (arquivo.bin)										   *
!*     SAVECTL - grava os descritores do "arquivo.bin" produzido por SAVEBIN   *
!*                                                                             *
!*     SAVEBIN_LINEAR - Grava campos de modelo em grade regular (arquivo.bin)  *
!*     SAVECTL_LINEAR - Grava os descritores do arquivo,bin produzido com      * 
!*                      SAVEBIN_LINEAR                                         *
!*                                                                             *
!*  ESTRUTURA DE DADOS PUBLICAS                                                * 
!*                                                                             *      
!*      XYZGRAD - ESTRUTURA DE CAMPO DE ALTITUDE COM MULTIPLOS NIVEIS          *
!*     stidtype - ESTRUTURA PARA DADOS DE ESTACAO METEOROLOGICA                *
!*                                                                             *
!*                                                                             *   
!*******************************************************************************
!* DEPENDENCIAS                                                                * 
!* MODULOS DATELIB E STRINGFLIB                                                * 
!*******************************************************************************
! HISTORICO 
!  DEZ 2002 -SHSF-   Versao Inicial ( Prototipo )
!  MAR 2006 -SHSF-   Inclusao/revisao das rotinas para ler e gravar em formato 
!                    de ponto de grades (Load/save)(bin/ctl)_linear. 
!                    Adaptacao de funcoes de interpolacao (1998 - MAER-CTA-IAR)

MODULE MGRADS


 USE DATELIB
 USE stringflib 
 implicit none

 PRIVATE



type stidtype
   character(len=8):: cod
   real::lat
   real::lon
   integer::nlev
end type  



type xyzgrad
   	real,pointer::var(:,:,:,:)
	character(len=4),POINTER::code(:)
	character(len=255),POINTER::name(:)
	integer,pointer::varlevs(:)        ! Numero de niveis de cada variavel
	integer::imax
	integer::jmax
	integer::kmax
	integer::nvars
	real,pointer::lev(:)
	real,pointer::lat(:)
	real,pointer::lon(:)
end type




	 
 PUBLIC stidtype
 PUBLIC XYZGRAD
 PUBLIC SAVECTL
 PUBLIC SAVECTL2
 PUBLIC SAVEBIN
 PUBLIC SAVEBIN_LINEAR
 PUBLIC SAVECTL_LINEAR
 PUBLIC LOADCTL_LINEAR
 PUBLIC LOADBIN_LINEAR
 PUBLIC INTERP2D
 PUBLIC INTERP3D

 logical::xrev ! Indica se o eixo x esta em ordem reversa
 logical::yrev ! Indica se o eixo y esta em ordem reversa   

 integer,parameter::recsize = 4 ! Tamanho da palavra ( 4 bytes ou 1  Word ) conforme compilador 

 
 CONTAINS
 

!==============================================================================!
! SaveCTL |                                                               |SHSF!
!==============================================================================!

 SUBROUTINE SAVECTL(un,filename,cdate,codes,desc,nvars)
      integer,intent(in):: un
      character(len=*),intent(in)::filename
      real*8,intent(in):: cdate
      character(len=*),dimension(:),intent(in)::codes ! Codigo das variaveis
      character(len=*),dimension(:),intent(in)::desc  ! Descricao das variaveis
      integer,intent(in)::nvars                       ! Numero das Variaveis

 
 !{ Variaveis locais 
	Character(len=255)::ctlname,filename2
	integer ::i,uni
        character(len=255)::gradsdate
	character(len=50),dimension(50)::substring
	integer ::nelements
!}

!{ Iniciar variaveis Locais
	uni=un
	i=0
        ctlname=trim(filename)//'.ctl'
!}	

!{ Cortar os diretorios do nome do arquivo
       call sep_substrings(filename,'/',substring,nelements)
       filename2=substring(nelements)
!} 

!{ Gravar arquivo descritor (CTL)
	open(uni,file=ctlname,status='unknown')
       
	write(uni,1001)trim(filename2)    
   1001 format('DSET ^',a,'.bin')
          
	    write(uni,1002)    
   1002 format('DTYPE station')
   
        write(uni,1003)trim(filename2)    
   1003 format('STNMAP ^',a,'.map')
    
        write(uni,1004) null
   1004 format('UNDEF ',f12.1)
   
        write(uni,1005)
   1005 format('TITLE INPE-CPTEC')

        gradsdate=grdate(cdate)
        write(uni,1006) trim(gradsdate)
   1006 format('TDEF 1 linear ',a,' 3hr')
   
        write(uni,1007)nvars
   1007 format('VARS ',i3)

        do i=1,nvars
          write(uni,'(a4,a5,1x,a30)')codes(i)," 0 99 ",desc(i)
        end do

        write(uni,'(a)')"ENDVARS"
!}
   close(uni)
 
 end subroutine  SAVECTL
 !==============================================================================!
! SaveCTL2 |                                                               |SHSF!
!==============================================================================!

 SUBROUTINE SAVECTL2(un,filename,cdate,codes,desc,nvars,levs,nlevs)
      integer,intent(in):: un
      character(len=*),intent(in)::filename
      real*8,intent(in):: cdate
      character(len=*),dimension(:),intent(in)::codes ! Codigo das variaveis
      character(len=*),dimension(:),intent(in)::desc  ! Descricao das variaveis
      integer,intent(in)::nvars                       ! Numero das Variaveis
	  real,dimension(:),intent(in)::levs              ! Niveis isobaricos
	  integer,intent(in)::nlevs                       ! Numero de niveis

 
 !{ Variaveis locais 
	Character(len=255)::ctlname,filename2
	integer ::i,uni
        character(len=255)::gradsdate
	character(len=50),dimension(50)::substring
	integer ::nelements,z
	character(len=1024)::XX
!}

!{ Iniciar variaveis Locais
	uni=un
	i=0
        ctlname=trim(filename)//'.ctl'
!}	

!{ Cortar os diretorios do nome do arquivo
       call sep_substrings(filename,'/',substring,nelements)
       filename2=substring(nelements)
!} 

!{ Gravar arquivo descritor (CTL)
	open(uni,file=ctlname,status='unknown')
       
	write(uni,1001)trim(filename2)    
   1001 format('DSET ^',a,'.bin')
          
	    write(uni,1002)    
   1002 format('DTYPE station')
   
        write(uni,1003)trim(filename2)    
   1003 format('STNMAP ^',a,'.map')
    
        write(uni,1004) null
   1004 format('UNDEF ',f12.1)
   
        write(uni,1005)
   1005 format('TITLE INPE-CPTEC')
        
		XX='ZDEF    '//trim(strs(nlevs))//' LEVELS'
		do z=1,nlevs
		 XX=trim(XX//' '//strs(levs(z)))
		end do
		write(uni,'(a)')trim(XX)


        gradsdate=grdate(cdate)
        write(uni,1006) trim(gradsdate)
   1006 format('TDEF 1 linear ',a,' 3hr')
   
        write(uni,1007)nvars
   1007 format('VARS ',i3)

        do i=1,nvars
          write(uni,'(a4,1x,i2,a3,1x,a30)')codes(i),nlevs," 99 ",desc(i)
        end do

        write(uni,'(a)')"ENDVARS"
!}
   close(uni)
 
 end subroutine  SAVECTL2

!==============================================================================!
! SaveBIN | Grava dados na forma de ponto de estacao                      |SHSF!
!==============================================================================!
 SUBROUTINE SAVEBIN(un,filename,obs,nobs,STID,nvars)
 
	
	integer,intent(in)::un
        character(len=*),intent(in)::filename
        type(stidtype),dimension(:),intent(inout)::STID  
        real,dimension(:,:),intent(INOUT)::obs
        integer,intent(in)::nobs
        integer,intent(in)::nvars
        integer ::irec

!{ Variaveis Locais


	Integer:: NFLAG,NLEV,i,uni,j,ub
	REAL::TIM,rlat,rlon
	character(len=255)::outfile
        character(len=8)::C
        character(len=4)::C1,C2
!}


!{ Iniciando Variaveis Locais 
	NFLAG = 1
	NLEV = 1
	TIM=0.0
        uni=un
        outfile=trim(filename)//".bin"
!}      
  
!{ Verificar se o tamanho da matriz/ nobs esta correto
      ub=ubound(obs,1)
      if (ub<nobs) then 
        Print *,"Erro em mgrads:"
        print *," nobs ou obs com tamanho incorretos"
        print *," nobs =",nobs
        print *," tamanho(nobs)=",ub
        stop
      end if

!}
!{ Grava arquivo Binario   

  
	open(uni,file=outfile,STATUS='unknown',FORM='UNFORMATTED',access='DIRECT',recl=recsize) 
	irec=0
	print *,nobs
		Do i=1,nobs
		C=STID(i)%COD
		C1=C(1:4)
		C2=C(5:8)
		irec=irec+1;WRITE (uni,rec=irec)C1
		irec=irec+1;WRITE (uni,rec=irec)C2
		irec=irec+1;WRITE (uni,rec=irec)STID(i)%Lat
		irec=irec+1;WRITE (uni,rec=irec)STID(i)%LON
		irec=irec+1;WRITE (uni,rec=irec)TIM
		irec=irec+1;WRITE (uni,rec=irec)NLEV
		irec=irec+1;WRITE (uni,rec=irec)NFLAG
		DO J=1,NVARS
			irec=irec+1
			WRITE (uni,rec=irec) OBS(I,J)
		END DO
	End do	

	rlat=0.0
	rlon=0.0
	TIM=0.0
	NLEV=0
	NFLAG=1
	irec=irec+1;WRITE (uni,rec=irec) C1
        irec=irec+1;WRITE (uni,rec=irec) C2
        irec=irec+1;WRITE (uni,rec=irec) rLat
        irec=irec+1;WRITE (uni,rec=irec) rLon
        irec=irec+1;WRITE (uni,rec=irec) TIM
        irec=irec+1;WRITE (uni,rec=irec) NLEV
        irec=irec+1;WRITE (uni,rec=irec) NFLAG
        close(uni)
        print *,trim(outfile)," ",irec*4," Bytes Gravados"

!}
	
END SUBROUTINE SAVEBIN





!==============================================================================!
! SaveBin_linear | produz campos no formato do grads                      |SHSF!
!==============================================================================!
!     Valores missing = -9.99e33 (Esto precisa ser revisto. Pode estar errado!)! 
! Esta subrotina grava campos de modelos no formato binario do grads           !            
! Os campos devem ser fornecidos em duas estruturas de dados:                  !  
!  a) Com os campos de superficie (xy)                                         !
!  b) Com os campos de altitude (xyz)                                          !
!                                                                              !  
!  O programa principal deve conter as declaracoes:                            !
!                                                                              !
!  use MGRADS                                                                  !
!  use datelib                                                                 !
!  type(xygrad):: xy         != Dados de superficie (x,y,nvar)                 !
!  type(xyzgrad) :: xyz   	 != Dados de altitude (x,y,z,nvars)                !
!                                                                              !
!                                                                              !
!       allocate(xy%code(NVar_Sup),STAT=ERR)                                   !
!		allocate(xy%name(NVar_Sup))                                            !
!		allocate(xyz%var(imax,jmax,kmax,nvar_alt))                             ! 
!		allocate(xyz%code(nvar_alt))                                           !
!		allocate(xyz%name(nvar_alt))                                           ! 
!		allocate(xyz%lev(nlevs))                                               !  
!                                                                              !
!   Onde Nvar_SUp  e o numero de campos (variaveis) de superficie              ! 
!        NVar_alt  e o numero de campos (variaveis) de altitude                !    								                                           !
!                                                                              ! 
!   Forneca em %code um codigo de 4 letras para cada variavel                  ! 
!   Forneca em %name o nome completo da variavel                               !
!   Forneca em %nlevs o  numero de niveis isobaricos dos campos de altitude    !
!   imax = numero de cordenadas x (longitudes)                                 !
!   jmax = numero de coordenadas y (latitudes)                                 !
!==============================================================================!
 subroutine savebin_linear(binname,xyz)

      
   !{  vaiaveis da interface 
   
    character(len=*),intent(in):: binname	! Nome do arquivo (sem extensao)
	type(xyzgrad)::xyz                      ! Dados de altitude (x,y,z,var)
	  
	!}
	  
	  
	  
	 !{- vari�eis locais	!
	  


	  integer ::x,y,z,v
	  integer*4 ::irec2
	  integer::nvars,i
	 
	  !}
	  !{ Inciando Vari�eis !
	  
      irec2=0;x=0;y=0;z=0;v=0;i=0;nvars=0
      
	 
      open(8,file=binname,status='unknown',FORM='UNFORMATTED',access='DIRECT',recl=recsize)

    

      	  
      !{ Salvando dados de superficie 
      
	  nvars=0 
	  

      write(*,*)'Saving...',trim(binname)
	  write(*,*)'imax,jmax,kmax =',xyz%imax,xyz%jmax,xyz%kmax
	  write(*,*)'nvar xy,xyz =',xyz%nvars
	  
	  
	      
	do v=1,xyz%nvars
		if (xyz%varlevs(v)>0) then 
		do z=1,xyz%varlevs(v)
			do y=1,xyz%jmax
				do x=1,xyz%imax
					irec2=irec2+1
					write(8,rec=irec2)xyz%var(x,y,z,v)		
				end do        
			end do             
		end do                
		else
		  z=1
			do y=1,xyz%jmax
				do x=1,xyz%imax
					irec2=irec2+1
					write(8,rec=irec2)xyz%var(x,y,z,v)		
				end do        
			end do             
		 endif	    

	end do                     
     
	print *,irec2*4,'Bytes'

	close(8)		
	end subroutine 
	 






!==============================================================================!
! SaveCtl_linear | produz descritor p/ SaveBin_linear                     |SHSF!
!==============================================================================!
	 

 SUBROUTINE SAVECTL_LINEAR(ctlname,binname,xyz,xmin,xmax,ymin,ymax,jdate)

 !{ Variaveis de interface
	character(len=*),intent(in)::ctlname
	character(len=*),intent(in)::binname
	type(xyzgrad),intent(in)::xyz
	
	real,intent(in)::xmin ! Longitude minima
	real,intent(in)::ymin ! 
	real,intent(in)::xmax ! 
	real,intent(in)::ymax ! 
	real*8,intent(in)::jdate ! Data em dias e fra�es de dias a partir do ano zero (vide datelib)
 !}

 !{ Variaveis locais

	character(len=255),dimension(10)::lines
	integer*2::z,i,v
	real :: stpx
	real :: stpy 
!}

	IF (XYZ%IMAX>0) stpx=(xmax-xmin)/real(xyZ%imax)
	IF (XYZ%JMAX>0) stpy=(ymax-ymin)/real(xyZ%jmax)
	  close(8)		
      open (8,file=(ctlname),status='replace')
    
        lines(1)='DSET    ^'//trim(binname)
		lines(2)='TITLE   Sample Data Set           '
        lines(3)='UNDEF   -9.99E33                  '
        lines(4)='XDEF    '//trim(strs(xyz%imax))//' LINEAR '//trim(strs(xmin))//' '//trim(strs(stpx))
        lines(5)='YDEF    '//trim(strs(xyz%jmax))//' LINEAR '//trim(strs(ymin))//' '//trim(strs(stpy))
        lines(6)='ZDEF    '//trim(strs(xyz%kmax))//' LEVELS'
		do z=1,xyz%kmax
		 lines(6)=trim(lines(6))//' '//trim(strS(xyz%lev(z)))
		end do

        lines(7)='TDEF    1 LINEAR '//grdate(jdate)
        lines(8)='VARS   '//trim(strS(xyz%nvars))

	do i=1,8
	    write(8,'(a)')trim(lines(i))
	end do

	do v=1,xyz%nvars
	    write(8,881)xyz%code(v),xyz%varlevs(v),trim(xyz%name(v))
 881    format(a4,i3,' 0 ',a)

	end do
   

        write(8,'(a)')'ENDVARS                            '
 
        close(8)	 


end subroutine







	 
!==============================================================================!
! LoadCtllinear | produz descritor p/ SaveBin_linear                      |SHSF!
!==============================================================================!
 SUBROUTINE LOADCTL_LINEAR(ctlname,binname,xyz)

 !{ Variaveis de interface
	character(len=*),intent(in)::ctlname
	character(len=*),intent(out)::binname
	type(xyzgrad),intent(out)::xyz


 !}

!{ Variaveis Locais
	character(len=255)::linha		     ! Uma linha de texto
	character(len=255),dimension(200)::w	 ! Palavras de uma linha de texto
	integer::nw                        ! Numero de palavras em w
	integer*2::i,II,jj
	integer::nlevs
	integer::err
	REAL::AUX


        YREV=.false.
	XREV=.false.

	  close(8)
      open (8,file=(ctlname),status='old')

  10	   read(8,'(a)',end=999) linha


          IF (index(ucases(linha),"YREV")>0) then
	       YREV=.true.
	       print *,"YREV=" ,YREV
	   end if
	  if (index(ucases(linha),"XREV")>0)  then
	      XREV=.true.
	      print *,"XREV=" ,XREV
	   end if

    	  !{ Separa nome do arquivo
	       if (index(linha,"DSET")>0)then
				call sep_words(linha,w,nw)
				binname=w(nw)
			end if
		  !}
		 !{ XDEF Definicao do eixo X
	       if (index(linha,"XDEF")>0)then
				call sep_words(linha,w,nw)
				XYZ%imax=val(w(2))
				if (associated(xyz%lon)) deallocate(xyz%lon)
				allocate(xyz%lon(xyz%Imax),STAT=ERR)
				if (err>0) stop
				IF (INDEX(W(3),"LEVELS")>0) THEN
					CALL GETLEVELS(W,NW,XYZ%imax,XYZ%LON)
				ELSE
					AUX=VAL(W(5))
					XYZ%LON(1)=VAL(W(4))
					DO II=2,XYZ%IMAX
						XYZ%LON(II)=XYZ%LON(II-1)+AUX
					END DO
				END IF

			
			end if
		  !}

		   !{ YDEF Definicao do eixo Y
	       if (index(linha,"YDEF")>0)then 
				call sep_words(linha,w,nw)
				XYZ%jmax=val(w(2))
				if (associated(xyz%lat)) deallocate(xyz%lat)
				allocate(xyz%lat(xyz%jmax),STAT=ERR)
				if (err>0) stop
				IF (INDEX(W(3),"LEVELS")>0) then 
					CALL GETLEVELS(W,NW,XYZ%jmax,XYZ%LAT)
				else
					AUX=VAL(W(5))
					XYZ%LAT(1)=VAL(W(4))
					DO II=2,XYZ%IMAX
						XYZ%LAT(II)=XYZ%LAT(II-1)+AUX
					END DO
				END IF
	
				
			end if
		  !}

		    !{ ZDEF Definicao do exio Z
	       if (index(linha,"ZDEF")>0)then 
				call sep_words(linha,w,nw)
				XYZ%kmax=val(w(2))

				!{ Obtendo niveis isobaricos 
				if (associated(xyz%lev)) deallocate(xyz%lev)
				allocate(xyz%lev(xyz%kmax),STAT=ERR)
				if (err>0) stop
				
				IF (INDEX(W(3),"LEVELS")>0) CALL GETLEVELS(W,NW,XYZ%Kmax,XYZ%LEV)

				
					
					 
			
					
		      !}

	
			end if
		  !}
		 !{ Separa nome do arquivo 
	       if (index(linha,"VARS")>0)then 
				
				call sep_words(linha,w,nw)
				xyz%nvars=val(w(2))
				if (associated(xyz%varlevs)) deallocate(xyz%varlevs)
				allocate(xyz%varlevs(xyz%nvars),STAT=ERR)

				if (associated(xyz%code)) deallocate(xyz%code)
				if (associated(xyz%name)) deallocate(xyz%name)
				allocate(xyz%code(xyz%nvars),xyz%name(xyz%nvars),STAT=ERR)
				if (ERR>0) stop 

				do i=1,xyz%nvars
						read(8,'(a)') linha
						call sep_words(linha,w,nw)
						xyz%varlevs(i)=val(w(2))
						xyz%code(i)=w(1)
					end do
					goto 999
				end if
			
		  !}

	  !{ 
	

		 goto 10
   999  continue 
		!{Aloca as estruturas xy e xyz
		!allocate(xy%code(NVar_Sup),STAT=ERR)
		!allocate(xy%name(NVar_Sup))
		if (associated(xyz%var)) deallocate(xyz%var)
		allocate(xyz%var(xyz%imax,xyz%jmax,xyz%kmax,xyz%nvars),STAT=ERR)
		
		if (err>0) stop
		

      close(8)	 

	end subroutine








!==============================================================================!
! LoadBin_linear | Carrega campos no formato do grads                     |SHSF!
!==============================================================================!
!     Valores missing = -9.99e33                                               !
! Esta subrotina carrega campos de modelos no formato binario do grads         !
! Os campos sao fornecidos em duas estruturas de dados:                        !
!  a) Com os campos de superficie (xy)										   !
!  b) Com os campos de altitude (xyz)                                          !
!                                                                              !  
!  O programa principal deve conter as declaracoes                             !:
!                                                                              !
!  use MGRADS                                                                  !
!  use datelib                                                                 !
!  type(xygrad):: xy         != Dados de superficie (x,y,nvar)                 !
!  type(xyzgrad) :: xyz   	 != Dados de altitude (x,y,z,nvars)                !
!                                                                              !
!==============================================================================!
! HISTORICO                                                                    !
!      DEZ 2002 -  S�gio H.S. Ferreira - prototipo                           !
!      MAR 2006 -  Sergio H.S. Ferreira - Revisao e adaptacao                  !
!==============================================================================!
 subroutine loadbin_linear(binname,xyz,v)

      
   !{  vaiaveis da interface 
   
    character(len=*),intent(in):: binname	! Nome do arquivo (sem extensao)
	type(xyzgrad),intent(inout)::xyz		! Dados de altitude (x,y,z,var)
	integer,intent(in):: v                  ! Indice da variavel a ser lida (1 a xyz%nvars)
	           
	  
	!}
	  
	  
	  
	 !{- vari�eis locais	!
	  

          real::bf ! Byte de formatacao 
	  integer ::x,y,z
	  integer*4 ::irec2
	  integer::nvars,i,i1,i2,is,j1,j2,js,k1,k2
	 
	  !}
	  !{ Inciando Vari�eis !
	  
      irec2=0;x=0;y=0;z=0;i=0;nvars=0
      
	  if ((v<1).or.(v>xyz%nvars))  stop
	
!      open(8,file=binname,status='unknown',FORM='UNFORMATTED',access='DIRECT',recl=(xyz%imax*xyz%jmax)*recsize)
	open(8,file=binname,status='unknown',access='DIRECT',recl=(xyz%imax*xyz%jmax+2)*recsize)
    
    !{ Caso o arquivo grads tena sido escrito em ordem reversa em x ou y, inverte a ordem dos indices 
       if (xrev) then 
          i1=xyz%imax;i2=1;is=-1
       else
          i1=1;i2=xyz%imax;is=1
       end if	
         
       if (yrev) then 
          j1=xyz%jmax;j2=1;js=-1
       else
          j1=1;j2=xyz%jmax;js=1
       end if	  
      !}
      
	  nvars=0 
	  
	  
	  
       
	    
	!}
	!{ Loading  Up Air and surface DATA 

       
	
		IREC2=VARREG(XYZ,V)
		
		
		if (xyz%varlevs(v)>0) then 
			k1=1
			k2=xyz%varlevs(v)
		else
			k1=1
			k2=1
		end if
		 do z=k1,k2
			irec2=irec2+1
			read(8,rec=irec2)bf,((xyz%var(x,y,z,v),x=i1,i2,is),y=j1,j2,js),bf
	 	end do
	
     

	close(8)	
	
		write(*,77)xyz%code(v),xyz%var(i1,j1,k1,v),xyz%lon(i1),xyz%lat(j1),xyz%lev(k1)
		write(*,78)xyz%var(i2,j2,k2,v),xyz%lon(i2),xyz%lat(j2),xyz%lev(k2)
77		format (1x,a,"[",f9.2,"(",f9.4,",",f9.4,",",f5.0,")->",\)	
78		format (f9.2,"(",f9.4,",",f9.4,",",f5.0,")]")
	end subroutine 
	 



!==============================================================================!
! itpl | Funcao interpolacao linear de dois pontos                        |SHSF!
!==============================================================================!
!  Funcao Privativa	 REAL
!  
!  Funcao simples que obtem a  interpolcacao de dois pontos Y1(x1), Y2(x2) 
!  no ponto x, isto e: 
!
!  Y=ITPL(X1,Y1,X2,Y2,X)
!
!      ^
!      |
!     Y2        +P2
!      |
!      Y-----+P
!      |     |
!     y1 +P1 |
!      |     |
!      --x1--x--x2-->
!
!    X1 = abscissa do ponto  P1
!    Y1 = ordenada do ponto  P1 
!    X2 = abscissa do ponto  P2
!    Y2 = ordenada do ponto  P2
!    X  = Abscissa do ponto que se deseja interpolar 
!    Y 	= Ordenada interolada
!    
!     Notas:
!     a) se X1<X<X2  Y e um ponto interpolado entre Y1 e Y2 
!        caso contrario Y sera um ponto extrapolado 
!	  b) Se X1=X2 entao Y=Y1=Y2
!=============================================================================



	Function ITPL(X1, Y1, X2, Y2, X)
		real:: ITPL
		real,intent(in)::X1,Y1 ! Ordenada e abscissa do ponto 1
		real,intent(in)::X2,Y2 ! Ordenada e abscissa do ponto 2
		real,intent(in)::X     ! Abscissa do ponto a ser interpolado(ou extrapolado)

		real :: ddd
	
		ddd = x2 - x1
		If (ddd /= 0) Then
			ITPL = (Y2 * (X - x1) + Y1 * (x2 - X)) / (x2 - x1) !: Def. Fun. Interpolacao  linear
		Else
			ITPL = Y1
		End If

	End Function  ITPL
	

!==============================================================================!
! INTERP2D | INTERPOLACAO ESPACIAL EM 2 DIMENSOES                         |SHSF!
!==============================================================================!
!  FUNCAO PUBLICA                  
!  
!  Obtem um valor interpolado de uma grade XY em um ponto de observacao
!  (Lat,Lon)
!
!  
!
!  Y=interp2d(VAR,VLAT,VLON,LAT,LON)
!	   
!    VAR (i=1:imax,j=1:jmax) = Campo (Pontos de grade) de uma variavel
!    Vlat(j=1:jmax)          = Latitudes de cada ponto "j" de Var 
!    Vlon(i=1:imax)          = Longitude de cada ponto "i" d Var(i,j)
!    Lat, Lon = Latitude e longitude para interpolacao 
!=============================================================================

	FUNCTION INTERP2D(Var,VLON,VLAT,lON,lAT);real::interp2d
	
	!{ Variaveis da interface
		real,dimension(:,:),intent(in)::Var ! Campo 2D do modelo (i,j) 
		real,dimension(:),intent(in):: Vlat ! Vetor de latitudes do campo Var Vlat(j)
		real,dimension(:),intent(in):: Vlon ! Vetor de longitudes do campo Var (Vlon(i))
		real,intent(in)::lat ! Latitude do ponto a ser interpolado
		real,intent(in)::lon ! Longitude do ponto a ser interpolado
	!}
	!{ Variaveis locais 
		integer::i,j,i1,j1,i2,j2
		integer::imax
		integer::jmax
		real::VX1,VX2 ! Variaveis auxiliares   
		
	 !}

		i=0;j=0
		i1=0;i2=0
		j1=0;j2=0
		INTERP2D=-9.9E33
		imax=ubound(Var,1)
		jmax=ubound(Var,2)
		
		
	!}

	!{ Verifica se os pontos nao correspondem aos estremos
	  if (vlat(1)>lat) then 
		j1=1
		j2=2
	  end if

	  if (vlat(jmax)<lat) then
		j1=jmax-1
		j2=jmax
	 end if

	 if (vlon(1)>lon) then
		i1=1
		i2=2
	 end if
	 if (vlon(imax)<lon) then
		i1=imax-1
		i2=imax
	end if

	!{ Localiza ponto de grade que contem as coordenada lat,lon 

		if (j1==0) then 
		do j=1,jmax-1
			  if ((vlat(j)<lat).and.(vlat(j+1)>=lat)) then 
				j1=j
				j2=j+1
				exit
			  end if
		end do
		if (j1==0) return 
		end if

		if (i1==0) then 
		do i=1,imax-1
		  

		  if ((vlon(i)<lon).and.(vlon(i+1)>=lon)) then 
				i1=i
				i2=i+1
				exit
		  	
			  end if
		end do
	 !}

		if (i1==0) return
		end if 

	  !{  Apos localizado a quadricula (i1,j1)-(i2,j2) que contem o ponto Lat,Lon, 
	  !   processa a interpolacao 
		if((i1*i2*j1*j2)>0) then 		
		  VX1=itpl(vlat(j1),var(i1,j1),vlat(j2),var(i1,j2),lat)
		  VX2=itpl(vlat(j1),var(i2,j1),vlat(j2),var(i2,j2),lat)			  
		  INTERP2D=itpl(vlon(i1),VX1,vlon(I2),VX2,LON)  
 		end if
 	  !}
	
	  end function


!==============================================================================!
! INTERP3D | INTERPOLACAO ESPACIAL EM 3 DIMENSOES                         |SHSF!
!==============================================================================!
!  FUNCAO PUBLICA                                                              ! 
!                                                                              ! 
!  Obtem um valor interpolado de uma grade XYz em um ponto de observacao       !
!  (Lat,Lon,lev)                                                               ! 
!                                                                              !
!                                                                              !
!                                                                              !
!  Y=interpdd(VAR,VLAT,VLON,LAT,LON)                                           ! 
!	                                                                           !
!    VAR (i=1:imax,j=1:jmax,k=1:kmax) = Campo (Pontos de grade) de uma variavel!
!    Vlat(j=1:jmax)          = Latitudes de cada ponto "j" de Var              !
!    Vlon(i=1:imax)          = Longitude de cada ponto "i" d Var(i,j)          !
!    vlev(k=1:kmax)          = Niveis isobaricos                               !
!    Lat, Lon = Latitude e longitude para interpolacao                         !    
!==============================================================================!

	FUNCTION INTERP3D(Var,VLON,VLAT,VLEV,lON,lAT,LEV);real::interp3d


	!{ Variaveis da interface
		real,dimension(:,:,:),intent(in)::Var ! Campo 2D do modelo (i,j) 
		real,dimension(:),intent(in):: Vlat ! Vetor de latitudes do campo Var Vlat(j) (graus)
		real,dimension(:),intent(in):: Vlon ! Vetor de longitudes do campo Var (Vlon(i)) (graus)
		real,dimension(:),intent(in)::Vlev ! Vetor de niveis isobaricos (hPa) 
		                                   ! (ordem decrescente de valores de pressao)
		real,intent(in)::lat ! Latitude do ponto a ser interpolado
		real,intent(in)::lon ! Longitude do ponto a ser interpolado
		real,intent(in)::lev ! Nivel isobarico (hPa)
	!}
	!{ Variaveis locais 
		integer::i,j,k,i1,j1,i2,j2,K1,K2
		integer::imax
		integer::jmax
		integer::kmax
		real::VX1,VX2,VZ1,VZ2 ! Variaveis auxiliares   
		
	 !}

		i=0;j=0;k=0
		i1=0;i2=0;k1=0
		j1=0;j2=0;k2=0
		INTERP3D=-9.9E33
		imax=ubound(Var,1)
		jmax=ubound(Var,2)
		kmax=ubound(var,3)

		
	!}

	!{ Verifica se os pontos nao correspondem aos extremos
	  if (vlat(1)>lat) then 
		j1=1
		j2=2
	  end if

	  if (vlat(jmax)<lat) then
		j1=jmax-1
		j2=jmax
	 end if

	 
	 if (vlon(1)>lon) then
		i1=1
		i2=2
	 end if
	 if (vlon(imax)<lon) then
		i1=imax-1
		i2=imax
	end if
	if (vlev(1)<vlev(2)) then
	   print *,"Erro! Niveis tem que ser informado em ordem descrescente" 
	   stop
	 end if

	 if (vlev(1)<lev) then
		k1=1
		k2=2
	 end if

	 if (vlev(kmax)>lev) then
		k1=kmax-1
		k2=kmax
	end if


	!{ Localiza ponto de grade que contem as coordenada lat,lon 

		if (j1==0) then 
		do j=1,jmax-1
			  if ((lat>vlat(j)).and.(lat<=vlat(j+1))) then 
				j1=j
				j2=j+1
				exit
			  end if
		end do
		if (j1==0) return 
		end if

		if (i1==0) then 
		do i=1,imax-1
		  

		  if ((lon>vlon(i)).and.(lon<=vlon(i+1))) then 
				i1=i
				i2=i+1
				exit
		  	
			  end if
		end do
	 !}

		if (i1==0) return
		end if 

	   if (k1==0) then
		do k=1,kmax-1
			if ((lev<=vlev(k)).and.(lev>vlev(k+1))) then 
				k1=k
				k2=k+1
				exit
			end if
		end do
		if(k1==0) return
		end if


	  !{  Apos localizado o cubo (i1,j1,K1)-(i2,j2,k2) que contem o ponto Lat,Lon,lev 
	  !   processa a interpolacao 
	
	 
	  
	  if((i1*j1*k1*i2*j2*k2)>0) then 	
	     !{ Obtem o valor interpolado da base do cubo
		VX1=itpl(vlat(j1),var(i1,j1,k1),vlat(j2),var(i1,j2,k1),lat)
		VX2=itpl(vlat(j1),var(i2,j1,k1),vlat(j2),var(i2,j2,k1),lat)			  
		VZ1=itpl(vlon(i1),VX1,vlon(I2),VX2,LON) 
		
		
	     !}

	     !{ Caso  seja o nivel isobarico exato  nao faz a interpolacao vertical
	     	if (vlev(k1)==lev) then 
			interp3d=VZ1
			
			RETURN
		END IF

	     !{ Obtem o valor interpolado da topo do cubo 
 		
		VX1=itpl(vlat(j1),var(i1,j1,k2),vlat(j2),var(i1,j2,k2),lat)
		VX2=itpl(vlat(j1),var(i2,j1,k2),vlat(j2),var(i2,j2,k2),lat)			  
		VZ2=itpl(vlon(i1),VX1,vlon(I2),VX2,LON) 


             !}

	     !{ Obtem o valor interoloado (interpolacao logaritima)

	     INTERP3D=itpl(log(vlev(k1)),VZ1,log(vlev(k2)),VZ2,log(lev))
           end if

 	  !}
	
	  end function



 


	  ! Esta subrotina e chamada exclusivamente por LOADCTL_LINEAR, PARA
	  ! fazer a leitura de XDEF, YDEF e ZDEF no caso destas nao serem 
	  ! definidas como "LINEAR" e sim como "LEVES"
	  !
	  !  W e NW corresponde ao resto de linha anterior e poderam conter 
	  !  alguns niveis caso NW > 3. Neste caso e feito o resto da leitura
	  !  destas variaveis. Caso ainda existam variaveis a serem lidas entao 
	  !  continua-se o processo 
	  !
	  !  
	  SUBROUTINE GETLEVELS(W,NW,JMAX,LVAR)
		character(len=*),dimension(:),intent(inout)::W
		integer,intent(inout)::nw
		integer,intent(in)::jmax ! Numero maximo de niveis 
		real,dimension(:),intent(inout)::LVAR

		integer:: jj,II
		CHARACTER(LEN=255)::LINHA

	
					JJ=0
					IF (NW>3) THEN
						DO II=4,NW
							JJ=JJ+1
							LVAR(JJ)=VAL(W(II))
						END DO
					end if

					DO WHILE (JJ<JMAX)
						read(8,'(a)')LINHA
						call sep_words(linha,w,nw)
						do ii=1,nw
							JJ=JJ+1
							LVAR(JJ)=val(w(ii))
						end do

					END DO
					 
	END SUBROUTINE

	 
	

	
!==============================================================================!
! VARREG | REGISTRO INICIAL DOS DADOS DE UMA VARIAVEL                     |SHSF!
!==============================================================================!
!  FUNCAO PRIVADA                                                              !
!                                                                              !
!  Calcula a posicao de uma variavel dentro de um arquivo grads, com base      !
!  nas informacoes lidas por LOADCTL_LINEAR ( ESTRUTURA xyz ) e indice         !
!  que identifica a variavel                                                   !
!                                                                              ! 
!                                                                              !
!==============================================================================! 
	function VARREG(XYZ,V);INTEGER*4::VARREG

	!{ Variaveis de interface
		TYPE(XYZGRAD),INTENT(IN)::XYZ  ! Estrutura de dados lidas por LOADCTL_LINEAR
		INTEGER,INTENT(IN)::V		   ! Indice da variavel 
	!}
	!{ Variaveis locais	 
		INTEGER*4::S
		integer :: II
	!}


	 IF ( V>1 )THEN
	    S=0 
		DO II=1,V-1
			IF (XYZ%VARLEVS(II)==0) THEN
				S=S+1
			ELSE
				S=S+XYZ%VARLEVS(II)
			END IF
		END DO
		    
			VARREG=S
	ELSE
			VARREG=0
	END IF

	END FUNCTION
				



	
END MODULE