!#######################################################################!
!  BRAMS/BAM - CPTEC/INPE - WRF-Chem - FIM-Chem emission models         ! 
!  version 1.8.3: Sep 2017                                              !
!  Coded by Saulo Freitas and Karla Longo                               !
!  brams support: brams_help@cptec.inpe.br - http://brams.cptec.inpe.br !
!#######################################################################!

module edgar_emissions
!---------------------------------------------------------------------------
integer, parameter :: maxnspecies = 200 , nspecies = 14

integer, parameter ::                 &
       CO          = 1  	      &
    ,  NOX         = 2                &
    ,  CO2         = 3                &
    ,  CH4	   = 4                &
    ,  SO2         = 5  	      &
    ,  N2O         = 6		      &
    ,  SF6	   = 7		      &
    ,  NMVOC	   = 8                &
    ,  SO4         = 9                &
    ,  URBAN2	   = 10               & !PM25
    ,  URBAN3      = 11               & !PM10
    ,  BC	   = 12               & 
    ,  OC          = 13               & 
    ,  NH3         = 14                 

!---------------------------------------------------------------------------
character(LEN=20),dimension(nspecies),parameter :: spc_name= &
! '12345678901234567890'
(/                         &
  'CO                  '   &	
, 'NOX                 '   &
, 'CO2                 '   &
, 'CH4                 '   &
, 'SO2                 '   &
, 'N2O                 '   &
, 'SF6                 '   &
, 'NMVOC               '   &
, 'SO4                 '   &
, 'URBAN2              '   &
, 'URBAN3              '   &
, 'BC                  '   &
, 'OC                  '   &
, 'NH3                 '   &
/)

real, dimension(:,:,:,:),allocatable :: RAWsrc

!---------------------------------------------------------------------------

  type edgar_vars   
     real, pointer, dimension(:,:,:)  :: src
!-----------

  end type edgar_vars

  type (edgar_vars), allocatable :: edgar_g(:)

contains
  !---------------------------------------------------------------

  subroutine alloc_edgar(edgar,n1,n2,n3)

    implicit none

    type (edgar_vars),dimension(nspecies)  :: edgar
    integer,intent(in) :: n1,n2,n3
    integer ispc
    
    do ispc=1,nspecies
     allocate (edgar(ispc)%src(n1,n2,n3))
    enddo

    return
  end subroutine alloc_edgar

  !---------------------------------------------------------------

  subroutine nullify_edgar(edgar)

    implicit none

    type (edgar_vars),dimension(nspecies)  :: edgar
    integer ispc

    do ispc=1,nspecies
       if (associated(edgar(ispc)%src))    nullify (edgar(ispc)%src)
    enddo

    return
  end subroutine nullify_edgar

end module edgar_emissions

!---------------------------------------------------------------
!---------------------------------------------------------------
  subroutine mem_edgar(n1,n2,n3)
    use edgar_emissions
    implicit none
    integer i
    integer, intent(in) :: n1,n2,n3

    if(.not. allocated(edgar_g)) allocate(edgar_g(nspecies))
    !do i=1,nspecies
    ! if(associated(edgar_g(i)%src)) deallocate(edgar_g(i)%src)
    !enddo

    call nullify_edgar(edgar_g(:))      
    call alloc_edgar  (edgar_g(:),n1,n2,n3) 
  end subroutine mem_edgar

!---------------------------------------------------------------
!---------------------------------------------------------------
subroutine read_edgar_antro(ihour,rhour,iyear,imon,iday,ng,ngrids,n1,n2,n3,rlat&
                           ,rlon,rland,deltax,deltay&
                           ,xt,yt,xm,ym,plat,plon)

use edgar_emissions

use grid_dims_out, only:            &
grid_type,                          &
edgar_data_dir,                     & 
fim_data_dir,                       & 
user_data_dir,                      &
use_edgar,                          &
cites_mobile_urban_emissions_file 

use util_geometry, only: &
T_mobileSrc,             &
loadMobileSrcMap,        &
polyToModel

use mem_grid, only :grid_g

implicit none

integer :: nlon    
integer :: nlat   
  
integer :: nmonths 
real   :: ilatn   
real   :: ilonn   


integer, intent (in) :: iyear
integer, intent (in) :: imon
integer, intent (in) :: iday
integer, intent (in) :: ng
integer, intent (in) :: n1
integer, intent (in) :: n2
integer, intent (in) :: n3
integer, intent (in) :: ngrids
integer, intent (in) :: ihour
integer, intent (in) :: rhour
!integer :: flag_inventario_RJ    !MArcelo 130516
!integer, dimension(5) :: pontos_MRJ    ! Angel 03061616   !pontos de cada uno de los 5 tipos de vias.
   

real, intent (in), dimension(n1,n2) :: rlat
real, intent (in), dimension(n1,n2) :: rlon
real, intent (in), dimension(n1,n2) :: rland
!real :: CO_MRJ    !Angel 030616
!real :: NOX_MRJ    !Angel 030616
!real :: old_emis    !Angel 030616
!real :: CO_via     !Angel 0306516
!real :: NOX_via     !Angel 030616
!real :: NOx_original_asc     !Angel DC 0306616
!real, dimension(5) :: porcentage    !Angel 030616 %de emision x tipo de via
!real, allocatable, dimension(:,:,:)  :: antiga_emision_nox     !Angel DC 030616


real, intent (in) :: deltax
real, intent (in) :: deltay

real,intent(in) :: xt(n1)
real,intent(in) :: yt(n2)
real,intent(in) :: xm(n1) 
real,intent(in) :: ym(n2)

character(len=250) :: prefix
character(len=250) :: suffix
character(len=250) :: filename(nspecies)
character(len=180) :: dummy

real,dimension(:),allocatable :: longedgar
real,dimension(:),allocatable :: latedgar

real :: plat
real :: plon

integer :: imonx
integer :: i
integer :: j  
integer :: ispc
integer :: nc
integer :: iread
integer :: im
integer :: ilaea
integer :: jlaea
integer :: fstat
integer :: i1
integer :: i2
integer :: j1
integer :: j2
integer :: ic
integer :: jc
integer :: k
   
real :: rrlat
real :: rrlon
real :: dlon1
real :: dlon2
real :: dlat1
real :: dlat2
real :: TX(nspecies)
real :: TX4(nspecies)
real :: lat
real :: lon
real :: rdummy
!real :: src_dummy(nmonths)

real, allocatable, dimension(:,:)     :: grid_area

real, allocatable, dimension(:,:) :: NOXAlonso
real, allocatable, dimension(:,:) :: COAlonso




!-------------------------------------------------------

!flag_inventario_RJ = 1 !InformaÃ§Ã£o do inventÃ¡rio local por vias para os jogos olimpicos !Angel 030616

!-------------------------------------------------------






type(T_mobileSrc), pointer, dimension(:) :: iMobileSrc

INTEGER :: edgar_ver
REAL :: ires,jres

integer :: flag_cuba_update
integer :: grade_cuba
integer :: termos_cuba

!Cuba
real, allocatable, dimension(:,:) :: NOXMade
real, allocatable, dimension(:,:) :: COMade
real, allocatable, dimension(:,:) :: NMVOCMade

integer, parameter ::  nlon_cuba = 130, nlat_cuba= 50, n_ps_cuba=8
real :: EMISS_CUBA,EMISS_CUBA2, EMISS_CUBA_PS, EMISS_CUBA_PS2
real, parameter:: ilatn_cuba=.1,ilonn_cuba=.1
real :: area_cuba
real :: accum_emission
integer :: ipoint,jpoint, ipoint_ANT, jpoint_ANT
character*240 prefix_c,suffix_c,filename2(nspecies)
! Stu variables for HTAP-FIM option 3/18/15
        integer :: ifim_case=0
        integer, allocatable, dimension(:,:) :: htap_fimN

flag_cuba_update = 0
grade_cuba = 0
termos_cuba=0
   
!   print*,"not doing EDGAR = returning" ; return
   
   
    
   ! SELECT CASE(use_edgar)
   ! CASE(1)
   !    edgar_ver = 3
   ! CASE(2)   
   !    edgar_ver = 4
   ! CASE(3)   
   !    edgar_ver = 5
   ! END SELECT
    	
    !SELECT CASE(edgar_ver)	
    IF(use_edgar==1) THEN
    ! CASE(3)
        nmonths = 1  		
        prefix='edg_'
        suffix='_2000_total.txt'
        ires = 1.
        jres = 1.
	nlon=360
	nlat=180
	ilonn=1.
	ilatn=1.
	allocate(longedgar(nlon))
	allocate(latedgar(nlat))
	if( .not. allocated (RAWsrc)) THEN
		allocate(  RAWsrc(nlon,nlat,nmonths,nspecies) )
		RAWsrc = 0.0
	endif

	do j=1,nlon;  longEdgar(j)=-179.5 + (j-1)*ilonn; enddo
	do i=1,nlat;   latEdgar(i)= -89.5 + (i-1)*ilatn; enddo
	
	if(ng == 1) then	
		call edgar3(nlon,nlat,nspecies,nmonths,spc_name,RAWsrc,edgar_ver,prefix,suffix,edgar_data_dir,ires,jres)
	endif



    !CASE(4 : 5)
    ELSEIF(use_edgar == 2 .or. use_edgar == 3) then 
        nmonths = 12	
        prefix='EDGAR'
        suffix='_2005.h5'	
        ires = .1
        jres = .1
	ilonn=0.1
	ilatn=0.1
	nlon=3600
	nlat=1800
	allocate(longedgar(nlon))
	allocate(latedgar(nlat))
	if( .not. allocated (RAWsrc)) THEN
		allocate(  RAWsrc(nlon,nlat,nmonths,nspecies) )
		RAWsrc = 0.0
	endif

	do j=1,nlon;  longEdgar(j)=-179.95 + (j-1)*ilonn; enddo
	do i=1,nlat;   latEdgar(i)= -89.95 + (i-1)*ilatn; enddo	
	
	if(ng==1) then
	   call edgar4(nlon,nlat,nspecies,nmonths,spc_name,RAWsrc,prefix,suffix,edgar_data_dir)
	endif
        !do ispc=1,nspecies
	! do imonx=1,12
	!  print*,"edgar4",ispc,imonx,maxval(RAWsrc(:,:,imonx,ispc))
        ! enddo
	!enddo

	IF( use_edgar == 3) then

          nmonths = 12
          PREFIX="EDGAR-HTAP_"
          SUFFIX="_2010.h5"
          ires = .1
          jres = .1
          ilonn=0.1
          ilatn=0.1
          nlon=3600
          nlat=1800
	  !----
          ifim_case=0
	  select case(grid_type)
	  case('fim')
           ifim_case=1
           if( .not. allocated (htap_fimN)) then
                  allocate(  htap_fimN(nlon,nlat) )
           endif
           open(11,FILE='htapLL_to_fimN.bin',FORM='UNFORMATTED')
           read(11)htap_fimN
           close(11)
           write(6,*)'htapLL_to_fimN.bin opened OK'
           CALL FLUSH(6)
	  !end select
           if( .not. allocated (grid_area)) then
          	  allocate(  grid_area(n1,n2) )
          	  grid_area = 0.0
           endif
	   call get_area_fim(grid_area,n1,n2) 
           accum_emission = 0.
           do i=1,n1
   		  accum_emission = accum_emission + grid_area(i,1)
           enddo
           write(6,*)'After get_area_fim,area sum=',accum_emission
           call flush(6)
	  end select
	  !----
          if( .not. allocated(longedgar)) allocate(longedgar(nlon))
          if( .not. allocated(latedgar )) allocate(latedgar(nlat))
          if( .not. allocated (RAWsrc)) then
          	  allocate(  RAWsrc(nlon,nlat,nmonths,nspecies) )
          	  RAWsrc = 0.0
          endif

          do j=1,nlon;  longEdgar(j)=-179.95 + (j-1)*ilonn; enddo
          do i=1,nlat;   latEdgar(i)= -89.95 + (i-1)*ilatn; enddo 

          write(6,*)'Before call to edgar_htap';CALL FLUSH(6)
          
	  if(ng==1) then
             call edgar_htap(nlon,nlat,nspecies,nmonths,spc_name,RAWsrc,prefix,suffix,edgar_data_dir)
          endif
          write(6,*)'After call to edgar_htap';CALL FLUSH(6)
	  
          do ispc=1,nspecies
	   do imonx=1,12
	     print*,"edgarhtap",ispc,imonx,maxval(RAWsrc(:,:,imonx,ispc))
           enddo
	  enddo

        ENDIF
    ELSE
    !CASE DEFAULT	
       STOP 'EDGAR version not allowed, check in: edgar_emissions.f90'
    !END SELECT  
    ENDIF
    
    
    !--- performs the interpolation to the model grid box	
    TX = 0.0
    do i=1,n1
   		do j=1,n2
    if(ifim_case.EQ.0)then 			
			call get_index1(i,j,nlon,nlat,n1,n2,rlat,rlon,longedgar,latedgar, &
		  			ilatn, ilonn ,dlat1,dlat2,dlon1,dlon2,i1,j1,i2,j2,ic,jc)

     			call interpol2(i,j,n1,n2,rlon,rlat,ic,jc,nlat,nlon,ilatn,ilonn,   &
	      			       imon,nmonths,nspecies,RAWsrc,tx(1:nspecies))
    else
        call grid_htap2fim(i,j,n1,n2,nlon,nlat,htap_fimN,imon,nmonths,nspecies,ilonn,ilatn, &
             latedgar,RAWsrc,grid_area,tx(1:nspecies))
    endif

 			! TX is the value interpolated to the model grid box.
			!-obs: the raw data is  the monthly mean and the units are in kg /(m^2 s)
     			do ispc = 1, nspecies 
     				edgar_g(ispc)%src(i,j,1)=TX(ispc)*86400. ! convert to kg/m2/day
     			enddo
                        !
			if(use_edgar == 1 .or. use_edgar == 2) then
			   !- special section for urban aerosols 
                           edgar_g(URBAN2)%src(i,j,1)=edgar_g(CO)%src(i,j,1)*(1./(28.*1e-3)) & ! mole[CO]/m2/day
                                                   *(2.95* 1.e-3) !=> kg[URBAN2]/m2/day 
                           edgar_g(URBAN3)%src(i,j,1)=edgar_g(CO)%src(i,j,1)*(1./(28.*1e-3)) & ! mole[CO]/m2/day
                                                   *(8.77* 1.e-3) !=> kg[PM10]/m2/day 
                          
			endif

   		enddo ! j =1,n2
    enddo ! i=1,n1




    if(trim(user_data_dir(1:len_trim(user_data_dir)) ) /= 'NONE' .and. &
        trim(user_data_dir(1:len_trim(user_data_dir)) ) /= 'none'  ) then
   		
		print*,'----------------------------------------------------------------------'
   		print*,'Using urban emissions updating for South America'
  	 	print*,'file=',trim(user_data_dir(1:len_trim(user_data_dir)))//'/'//trim(cites_mobile_urban_emissions_file)
		
		!------ calculate the grib box area, if necessary
          if( .not. allocated (grid_area)) then
   		allocate(grid_area(n1,n2))
          endif
		
		select case(grid_type)
		case('rams')
			call get_area_rams(grid_area,n1,n2,xt,yt,xm,ym)
		case('polar')
			call get_area_rams(grid_area,n1,n2,xt,yt,xm,ym)
		case('ll')
			call get_area_ll(grid_area,n1,n2)
		case('lambert')
			grid_area(:,:) = 1./(grid_g(ng)%dxt(:,:)*grid_g(ng)%dyt(:,:))
		case('mercator')
			grid_area(:,:) = 1./(grid_g(ng)%dxt(:,:)*grid_g(ng)%dyt(:,:))
		case('fim')
			call get_area_fim(grid_area,n1,n2) 
		case('gg')
			call get_area_gg(grid_area,n1,n2,rlat,rlon) 
		case default
			stop 'grid_type wrong'	
		end select
	
	 allocate(NOXAlonso(n1,n2))
   	 allocate(COAlonso(n1,n2))
   
   	 NOXAlonso = 0.0
   	 COAlonso  = 0.0
   
   	 call loadMobileSrcMap(trim(user_data_dir)//'/'//trim(cites_mobile_urban_emissions_file), iMobileSrc)
	
	 call polyToModel(iMobileSrc,ng,n1,n2,xt,yt,deltax,deltay,plat,plon,rlon,rlat,grid_type, COAlonso, NOXAlonso)
	
	 do j = 1, n2
   		do i = 1, n1
		
			COAlonso(i,j)  = ((COAlonso(i,j)*1e7)/365.) / grid_area(i,j)  ! kg/m2/dia
			NOXAlonso(i,j) = ((NOXAlonso(i,j)*1e7)/365.) / grid_area(i,j) ! kg/m2/dia
		
			if(NOXAlonso(i,j) .gt. 0)then
				do ispc=1,nspecies
					if (ispc.ne.NOX .and. ispc.ne.CO .and. edgar_g(NOX)%src(i,j,1) .ne. 0.) then
			        		edgar_g(ispc)%src(i,j,1) = (  NOXAlonso(i,j) / edgar_g(NOX)%src(i,j,1) )  * edgar_g(ispc)%src(i,j,1) ! nox ratio 
					end if 
				end do

				edgar_g(NOX)%src(i,j,1)  = NOXAlonso(i,j)
				edgar_g(CO)%src(i,j,1)   = COAlonso(i,j)
				!
				!srf- special section for urban aerosols 
                                if(use_edgar == 1 .or. use_edgar == 2) then
				  edgar_g(URBAN2)%src(i,j,1)=edgar_g(CO)%src(i,j,1)*(1./(28.*1e-3)) & ! mole[CO]/m2/day
                                                   *(2.95* 1.e-3) !=> kg[PM25]/m2/day 
                                  edgar_g(URBAN3)%src(i,j,1)=edgar_g(CO)%src(i,j,1)*(1./(28.*1e-3)) & ! mole[CO]/m2/day
                                                   *(8.77* 1.e-3) !=> kg[PM10]/m2/day 
                                endif

			endif
	
		enddo
   	 enddo
	
	 deallocate(NOXAlonso, COAlonso)
   	 deallocate(grid_area)

     endif ! if(trim(user_data_dir(1:len_trim(user_data_dir)) ) /= 'NONE'
	
     !- deallocate memory that we do not need anymore    
     if (ng==ngrids) deallocate (RAWsrc)
 
     if(grid_type == 'rams' .or. grid_type == 'polar') then ! only for 'rams' until rland is also defined for others grids 
    		do ispc=1,nspecies 
       			call apply_land_restriction(n1,n2,rland,edgar_g(ispc)%src(:,:,1))
    		enddo  
     endif
	
     !--------------------------------------Cuba Inventories--------------------------------

     if(flag_cuba_update == 1) then

		!------ calculate the grib box area, if necessary

         if(.NOT. allocated(grid_area)) allocate(grid_area(n1,n2))

		
		select case(grid_type)
		case('rams')
			call get_area_rams(grid_area,n1,n2,xt,yt,xm,ym)
		case('polar')
			call get_area_rams(grid_area,n1,n2,xt,yt,xm,ym)
		case('ll')
			call get_area_ll(grid_area,n1,n2)
		case('lambert')
			grid_area(:,:) = 1./(grid_g(ng)%dxt(:,:)*grid_g(ng)%dyt(:,:))
		case('mercator')
			grid_area(:,:) = 1./(grid_g(ng)%dxt(:,:)*grid_g(ng)%dyt(:,:))
		case('fim')
			call get_area_fim(grid_area,n1,n2) 
		case('gg')
			call get_area_gg(grid_area,n1,n2,rlat,rlon) 
		case default
			stop 'grid_type wrong'	
		end select
	
          if (grade_cuba==0) then


             !- name of dataset
             prefix_c='emis_cuba_'
             suffix_c='.csv'
           
             do ispc=1,nspecies

        	 if(trim(spc_name(ispc))=='URBAN2' .or. trim(spc_name(ispc))=='URBAN3') cycle

        	 nc=len_trim(spc_name(ispc))
        	 filename2(ispc)=trim(user_data_dir)//'/'//trim(prefix_c)//spc_name(ispc)(1:nc)//suffix_c
        	 print *,'opening   ', trim(filename2(ispc)),' - specie= ',trim(spc_name(ispc))
        	 open(12,file=trim(filename2(ispc)),status='old')

              !read  header

        	 do iread=1,nlon_cuba*nlat_cuba !maximum number of points for nlon,nlat
        	     EMISS_CUBA=0
        	     EMISS_CUBA2=0
        	    read(12,*,end=12) lon,lat,EMISS_CUBA,area_cuba
        	     rrlon = lon
        	     rrlat = lat
        	     call update_emissions_by_point(ng,n1,n2,xt,yt,deltax,deltay,plat,plon,rlon,rlat&
        				 ,rrlon,rrlat,ipoint,jpoint,grid_type)
        		    edgar_g(ispc)%src(ipoint,jpoint,1) = 0.
        		 if(ipoint /= -999 .and. jpoint /= -999) then
        		    EMISS_CUBA2 = ((EMISS_CUBA*1e6)/365)/(grid_area(ipoint,jpoint))
        		    edgar_g(ispc)%src(ipoint,jpoint,1) = EMISS_CUBA2	 
        		 endif
        	 enddo !of iread

        	 12 close(12)
             enddo
          else if (grade_cuba == 1) then
		
   	      allocate(NOXMade(n1,n2))
   	      allocate(COMade(n1,n2))
   	      allocate(NMVOCMade(n1,n2))
   	 
   	      NOXMade = 0.
   	      COMade  = 0.
   	      NMVOCMade  = 0.
   	 
   	      call loadMobileSrcMap(trim(user_data_dir)//'/Cuba_prov_Emissions.asc', iMobileSrc, 1 )
   	      
   	      call polyToModel(iMobileSrc,ng,n1,n2,xt,yt,deltax,deltay,plat,plon,rlon,rlat,grid_type, COMade, NOXMade, NMVOCMade)
   	      
   	      do j = 1, n2
   		      do i = 1, n1
   		      
   			      COMade(i,j)  = ((COMade(i,j)*1e7)/365.) / grid_area(i,j)  ! kg/m2/dia
   			      NOXMade(i,j) = ((NOXMade(i,j)*1e7)/365.) / grid_area(i,j) ! kg/m2/dia
   			      NMVOCMade(i,j) = ((NMVOCMade(i,j)*1e7)/365.) / grid_area(i,j) ! kg/m2/dia
   			      if(COMade(i,j) .gt. 0.)then
   				      edgar_g(NOX)%src(i,j,1)  = NOXMade(i,j)
   				      edgar_g(CO)%src(i,j,1)   = COMade(i,j)
   				      edgar_g(NMVOC)%src(i,j,1)   = NMVOCMade(i,j)
   			      end if
   				      !
   				      !srf- special section for urban aerosols 
				      if(use_edgar == 1 .or. use_edgar == 2) then
   				        edgar_g(URBAN2)%src(i,j,1)=edgar_g(CO)%src(i,j,1)*(1./(28.*1e-3)) & ! mole[CO]/m2/day
   							 *(2.95* 1.e-3) !=> kg[PM25]/m2/day 
   				        edgar_g(URBAN3)%src(i,j,1)=edgar_g(CO)%src(i,j,1)*(1./(28.*1e-3)) & ! mole[CO]/m2/day
   							 *(8.77* 1.e-3) !=> kg[PM10]/m2/day 
   				      endif
   	      
   		      end do
   	      end do
	      deallocate(NOXMade, COMade, NMVOCMade)	
          end if

   	   if (termos_cuba == 1) then
   	   !- name of dataset
   	   prefix_c='emis_cuba_ps_'
   	   suffix_c='.csv'
   	 

   	   do ispc=1,nspecies
   	   
   		 
   	 
   	     if(trim(spc_name(ispc))=='URBAN2' .or. trim(spc_name(ispc))=='URBAN3') cycle
   	 
   	     nc=len_trim(spc_name(ispc))
   	     filename2(ispc)=trim(user_data_dir)//'/'//trim(prefix_c)//spc_name(ispc)(1:nc)//suffix_c
   		 
   	     print *,'opening	', trim(filename2(ispc)),' - specie= ',trim(spc_name(ispc))
   	     open(13,file=trim(filename2(ispc)),status='old')
   	 
   	      !read  header
   	 
   	     do iread=1, n_ps_cuba !maximum number of points for nlon,nlat
   		     
   		 read(13,*,end=13) lon,lat,EMISS_CUBA
   		     rrlon = lon
   		     rrlat = lat
   		     call update_emissions_by_point(ng,n1,n2,xt,yt,deltax,deltay,plat,plon,rlon,rlat&
   					 ,rrlon,rrlat,ipoint,jpoint,grid_type)
   		     if(ipoint /= -999 .and. jpoint /= -999) then
   		      accum_emission = 0.
   		      if(ipoint_ANT==ipoint .and. jpoint_ANT==jpoint) accum_emission=edgar_g(ispc)%src(ipoint,jpoint,1)
   		      EMISS_CUBA2 = ((EMISS_CUBA*1e6)/365)/(grid_area(ipoint,jpoint))
   		      edgar_g(ispc)%src(ipoint,jpoint,1) = EMISS_CUBA2+accum_emission
            
                      ipoint_ANT=ipoint
                      jpoint_ANT=jpoint
                      endif
             enddo !of iread
              13 close(13)
           enddo
    end if

   deallocate(grid_area)
   deallocate(longedgar)
   deallocate(latedgar)

  endif

end subroutine read_edgar_antro
!------------------------------------------------------------------------------
subroutine edgar3(nlon,nlat,nspecies,nmonths,spc_name,RAWsrc,edgar_ver,prefix,suffix,edgar_data_dir,ires,jres)
implicit none

 integer, intent(in)::nlon,nlat,nspecies,nmonths,edgar_ver
 real, intent(in)::ires,jres
 character(len=*),intent(in)::spc_name(nspecies),prefix,suffix,edgar_data_dir
 real, intent(out)::RAWsrc(nlon,nlat,nmonths,nspecies)
 
 real::lat,lon,src_dummy(nmonths)
 character(len=240)::filename
 integer :: ispc,ilaea,jlaea,iread,im, fstat
 character *10 dummy




	do ispc=1,nspecies

		filename=trim(edgar_data_dir)//'/'//trim(prefix)//trim(spc_name(ispc))//suffix
		print *,'opening EDGAR emissions', trim(filename),' - specie= ',trim(spc_name(ispc))
		print*,'----------------------------------------------------------------------'
		if(trim(spc_name(ispc))=='URBAN2' .or. trim(spc_name(ispc))=='URBAN3' .or. &
		trim(spc_name(ispc))=='NH3') cycle
		
		open(11,file=trim(filename),status='old')
			!read  header
			read(11,*) dummy
			DO WHILE (dummy /= '#Values')
				read(11,*) dummy 
				!print*,dummy     
			enddo
       
		
			open(11,file=trim(filename),status='old')			
			DO iread=1,nlon*nlat !maximum number of points for nlon,nlat    
				read(11,*,iostat=fstat) lon,lat,(src_dummy(im),im=1,nmonths)
				if(fstat .ne. 0)exit
					!print*,lon,lat,(src_dummy(im),im=1,nmonths)
					!stop 33
					ilaea = nint(((lon - (-180.))/ires)+ 1 ) ! 2. = 1/resolution
					jlaea = nint(((lat - ( -90.))/jres)+ 1 )  ! 0.25 = 0.5*resolution
					ilaea= max(ilaea,1)

					if(ilaea < 1 .or. ilaea > nlon) stop 33
				
					if(jlaea < 1 .or. jlaea > nlat) stop 35
     				
					RAWsrc(ilaea,jlaea,:,ispc)=src_dummy(:)	
				
				
			END DO
			close(11)

		

	enddo

end subroutine edgar3
!------------------------------------------------------------------------------------
!> @author
!> GMAI
!> @version
!> 1.0
!> @date
!> 05/04/2012
!> @brief
!> DESCRIPTION
!> Reading function of the hd5 input file from EDGAR version 4.1 with attributes and &
!> organized by emission sources (stationary sources, mobile sources, aviation sources) to be read by the program PREP-CHEM-SRC.

!----------------------------------------------------------------------------------
subroutine edgar4(nlon,nlat,nspecies,nmonths,spc_name,RAWsrc,prefix,suffix,edgar_data_dir)
use hdf5
implicit none
 !variaveis de entrada e saída
 integer, intent(in):: nlon,nlat,nspecies,nmonths
 character(len=*), intent(in)::spc_name(nspecies)
 real,intent(out):: RAWsrc(nlon,nlat,nmonths,nspecies)
 character(len=*), intent(in)::prefix,suffix,edgar_data_dir
 
!variaveis do hdf5
 character(len=50) ::dsetname
 integer(hid_t):: file_id, dset_id, group_id, dtype_id
 integer(hsize_t) :: data_dims(2) 
 character(len=240)::filenameed
 integer :: error, printflags,ispc,ifontes
 real,dimension(:,:),allocatable::soma
 real,dimension(:,:,:),allocatable::somat
 integer(hid_t)::dcpl
 integer,dimension(1:1)::cd_values
 integer(size_t)::nelmts
 integer:: flags, filter_id, filter_info_both,filter_info
 integer, parameter:: fontes=3
 integer(size_t),parameter:: MaxChrlen=80
 character(len=MaxChrlen):: name
 integer, parameter :: STATIONARY_SOURCES = 1, &
	MOBILE_SOURCES = 2, &
	AVIATION_SOURCES = 3
 character(LEN=19),dimension(fontes), parameter :: &
	groupname=(/ &
	'STATIONARY_SOURCES', &
	'MOBILE_SOURCES    ', &
	'AVIATION_SOURCES  '/)
 integer :: imonx

  allocate(soma(nlon,nlat))
  allocate(somat(nlon,nlat,nmonths))

  !variavel responsavel por mostrar mensagens de erro do hdf5, valor 0 implica que sua funcao esta desativada
  printflags=0
	
  filenameed=trim(edgar_data_dir)//'/'//trim(prefix)//suffix
  print*, trim(filenameed)
  !inicializa interface do hdf5
  call h5open_f(error)
  print*, 'After h5open_f call'
    if(error==0)then
      !nao exibe mensagem de erro do hdf5 caso ocorra um
      call h5eset_auto_f(printflags, error)
      !abre arquivo com extensao hdf5
      call h5fopen_f(filenameed, H5F_ACC_RDONLY_F, file_id, error)
        if(error==0)then
          !nao exibe mensagem de erro do hdf5 caso ocorra um
          call h5eset_auto_f(printflags, error)
          !loop referentes as especies
          do ispc=1,nspecies
            soma=0.
            somat=0.
            dsetname=trim(spc_name(ispc))
            
            if(trim(spc_name(ispc))=='URBAN2' .or. &
               trim(spc_name(ispc))=='URBAN3' .or. trim(spc_name(ispc))=='SO4') cycle
            
            do ifontes=1, fontes
             
              print*,'----------------------------------------------------------------------'
	      print *,'opening EDGAR emissions', trim(filenameed),' - group=',trim(groupname(ifontes)),&
                ' - specie= ',trim(spc_name(ispc))
             
	      !abre o grupo fontes fixas dentro arquivo hdf5
              call h5gopen_f(file_id, groupname(ifontes), group_id, error)
                if(error==0) then
                  !nao exibe mensagem de erro do hdf5 caso ocorra um
                  call h5eset_auto_f(printflags, error)
                  !abre o dsetname(dsetname=especie)
                  call h5dopen_f(group_id, dsetname, dset_id, error)
                    if(error==0) then
                    !nao exibe mensagem de erro do hdf5 caso ocorra um
                      call h5eset_auto_f(printflags, error)
                      !captura o valor do dset
                      call h5dget_create_plist_f(dset_id,dcpl,error)
                        nelmts = 1
                        call H5Pget_filter_f(dcpl, 0, flags, nelmts,&
                        cd_values, MaxChrLen, name, filter_id, error)
                        if(error==0) then
                          !nao exibe mensagem de erro do hdf5 caso ocorra um
                          call h5eset_auto_f(printflags, error)
                          !leitura da soma de fontes fixas
                          call h5dread_f(dset_id, H5T_IEEE_F32LE, soma,&
                          data_dims, error)
                            if(error==0) then
                              !nao exibe mensagem de erro do hdf5 caso ocorra um
                              call h5eset_auto_f(printflags, error)
                
                              print *,dsetname,' ',maxval (soma),minval (soma)
                            
                              somat(:,:,1)=somat(:,:,1)+soma(:,:)
                              
                            endif
                        endif
                      !fecha o dcpl
                      call h5pclose_f(dcpl,error)
                    endif
                  !fecha o dset
                  call h5dclose_f(dset_id,error)
                endif
      
              !fecha o grupo fontes fixas
              call h5gclose_f(group_id, error)
            enddo
	    do imonx = 1,nmonths
              RAWsrc(:,:,imonx,ispc)=somat(:,:,1)
	    enddo
          enddo
          
          deallocate(soma)
          deallocate(somat)
        else
          print*,'EDGAR4 file not found: ', trim(filenameed)
          print*,'check the path of this HDF5 file'
          stop "EDGAR4 routine 1"
        endif
      call h5fclose_f(file_id, error)
    else 
      print*,'HDF5 library files not found'
      stop "EDGAR4 routine 2"
    endif
    !fecha a inicialização do hdf5
  call h5close_f(error)



end subroutine edgar4

!---------------------------------------------------------------
subroutine  edgar_htap(nx,ny,nspecies,nmonths,spc_name,RAWsrc,prefix,suffix,edgar_data_dir)


 use hdf5
 implicit none
 integer, intent(in) :: nx,ny,nmonths,nspecies
 real,intent(inout):: RAWsrc(nx,ny,nmonths,nspecies)
 character(len=*), intent(in)::prefix,suffix,edgar_data_dir
 character(len=*), intent(in)::spc_name(nspecies)

 !integer, parameter :: NX=3600
 !integer, parameter :: NY=1800
 !integer, parameter :: NZ=1
 !integer, parameter :: NT=12
 integer, parameter :: NSETOR=5
 !integer, parameter :: NSPECIES=9
 integer :: i
 integer :: j
 !integer :: imonth
 integer :: isetor
 integer :: ispc
 real, dimension(:,:,:), allocatable :: src_dummy
 character(len=11), dimension(NSETOR) :: setor=(/&
	'AGRICULTURE',&
	'ENERGY     ',&
	'INDUSTRY   ',&
	'RESIDENTIAL',&
	'TRANSPORT  '/) 
! character(len=5), dimension(NSPECIES) :: specie=(/&
!	'BC',&
!	'CO',&
!	'NH3',&
!	'NMVOC',&
!	'NOx',&
!	'OC',&
!	'PM10',&
!	'PM2.5',&
!	'SO2'/)
! character(len=25), parameter :: PREFIX="EDGAR-HTAP_"
! character(len=25), parameter :: SUFIX="_2010.h5"
! character(len=255), parameter :: DIR="/scratchin/grupos/catt-brams/home/valter.oliveira/PROJETOS/MATRIX/src/dados/"


 character(len=255) :: filename
 logical :: file_exists
 !variaveis do hdf5
 integer(HID_T) :: file_id
 integer(HID_T) :: dset_id 
 integer(HSIZE_T) :: data_dims(4)
 integer :: hdferr
 character(len=25) :: dsetname


   !- if a specie is included in this HTAP inventory, it will be written over the previously
   !- quantity defined by EDGAR 4 inventory
   do ispc=1, NSPECIES !loop de especies
   
     !- species not included in this inventory
     if(trim(spc_name(ispc))=='CO2' .or. &
     	trim(spc_name(ispc))=='CH4' .or. &
     	trim(spc_name(ispc))=='SO4' .or. &
     	trim(spc_name(ispc))=='N2O' .or. &
     	trim(spc_name(ispc))=='SF6'  )  then
        cycle
     else
     	RAWsrc(:,:,:,ispc) = 0.0
     endif   
   enddo   
   !
 

   allocate(src_dummy(NX,NY,nmonths))
   call h5open_f(hdferr) !inicializa o hdf5

   do isetor=1, NSETOR !loop de setores
      
      filename=trim(edgar_data_dir)//"/"//trim(PREFIX)//trim(setor(isetor))//trim(suffix)
      !print*,"filename=",filename;pause 333
      inquire(file=trim(filename), exist=file_exists)
      if(file_exists) then
         call h5fopen_f(trim(filename),H5F_ACC_RDONLY_F, file_id, hdferr)  !abre o arquivo
            
            do ispc=1, NSPECIES !loop de especies
               !dsetname=trim(specie(ispc))
               dsetname=trim(spc_name(ispc))
               
	       !- species not included in this inventory => cycle
	       if(trim(spc_name(ispc))=='CO2' .or. &
	          trim(spc_name(ispc))=='CH4' .or. &
	          trim(spc_name(ispc))=='SO4' .or. &
	          trim(spc_name(ispc))=='N2O' .or. &
	          trim(spc_name(ispc))=='SF6'  ) cycle
	       
	       !- species included in this inventory with different name
	       if(trim(spc_name(ispc))=='URBAN2')dsetname="PM2.5"
	       if(trim(spc_name(ispc))=='URBAN3')dsetname="PM10"
               
	       !- agriculture inventory only has inventory for NH3 
	       if(trim(setor(isetor)) .eq. 'AGRICULTURE' .and. trim(spc_name(ispc)) .ne. 'NH3') cycle
                
               call h5dopen_f(file_id, trim(dsetname), dset_id, hdferr)  !abre o dataset

               call h5dread_f(dset_id, H5T_NATIVE_REAL, src_dummy, data_dims, hdferr)  !faz a leitura do dataset
               
               do i=1,nmonths
                  print*,trim(setor(isetor)),"   ",trim(dsetname),i, " MAX ",maxval(src_dummy(:,:,i))
                  print*,trim(setor(isetor)),"   ",trim(dsetname),i, " MIN ",minval(src_dummy(:,:,i))                    
                  !-accumulate over sectors
		  RAWsrc(:,:,i,ispc) = RAWsrc(:,:,i,ispc) + src_dummy(:,:,i)
	       enddo

               call h5dclose_f(dset_id,hdferr)  !fecha o dataset
                
            enddo
         call h5fclose_f(file_id, hdferr)  !fecha o arquivo

      else
         print*, "EDGAR_HTAP file not found: ",trim(filename)
	 stop 
      endif
   enddo
 call h5close_f(hdferr)
 deallocate(src_dummy)
end subroutine  edgar_htap
!---------------------------------------------------------------
!---------------------------------------------------------------
subroutine get_edgar_indentity(spc_name,ident)
!use chem1_list
use edgar_emissions, only :  edgar_nspecies=>nspecies&
                            ,edgar_spc_name=>spc_name&
, edgar_CO     => CO	  & ! 1
, edgar_NOX    => NOX	  & ! 2  
, edgar_CO2    => CO2	  &   
, edgar_CH4    => CH4	  &  
, edgar_SO2    => SO2	  & 
, edgar_N2O    => N2O	  &
, edgar_SF6    => SF6	  &
, edgar_NMVOC  => NMVOC   &
, edgar_SO4    => SO4     &
, edgar_URBAN2 => URBAN2  &
, edgar_URBAN3 => URBAN3  &
, edgar_BC     => BC      &
, edgar_OC     => OC      &
, edgar_NH3    => NH3       ! 14


implicit none
integer isp
character (len=*), intent(in)  :: spc_name
integer          , intent(out) :: ident

do isp = 1,edgar_nspecies
  ident=-1
  if(spc_name == edgar_spc_name(isp)) then
      print*,'==>edgar found for ',spc_name
      ident=isp
      return
   endif
enddo

!print*,'chem1-list specie ',trim(spc_name), ' does not match if any one of edgar'
!stop 444
end subroutine get_edgar_indentity
!---------------------------------------------------------------
        subroutine grid_htap2fim(i,j,n1,n2,nlon,nlat,htap_fimN,imon,nmonths,nspecies,ilonn,ilatn, &
             latedgar,RAWsrc,grid_area,tx)
implicit none
integer, dimension(nlon,nlat) :: htap_fimN
integer nlat,nlon,i,j,imon,nmonths,nspecies,n1,n2
real, dimension(nlon,nlat,nmonths,nspecies) :: RAWsrc
real, dimension(n1,n2) :: grid_area
real, dimension(nlat) :: latedgar
real tx(nspecies)
real ilonn,ilatn
!-local var
integer ncount,ii,jj
real dlonr,dlatr,d2r,rearth
   rearth=6371.67e3  ! earth radius in meters, need in same units as variable grid_area
   d2r = 4.0*ATAN(1.0)/180.0
   dlonr=d2r*(ilonn)*rearth
   dlatr=d2r*(ilatn)*rearth
ncount = 0
TX(:)  = 0.
do jj = 1,nlat
   do ii = 1,nlon
    if(htap_fimN(ii,jj).eq.i)then
     	    ncount = ncount + 1
! Convert from kg/m(2)/sec to kg/sec and sum over FIM grid = i
     	    TX(:) = TX(:) + cos(latedgar(jj)*d2r)*dlonr*dlatr*RAWsrc(ii,jj,imon,:)  
    endif
   enddo
enddo
if(ncount.eq.0) STOP 'NCOUNT=0 in grid_htap2fim'
! Convert back to kg/m(2)/sec - now normalized over FIM grid
TX(:) = TX(:) / grid_area(i,j)
end subroutine grid_htap2fim
!---------------------------------------------------------------
