!###############################################################################!
!  CCATT-BRAMS/MCGA-CPTEC/WRF emission model    CPTEC/INPE                      !
!  Version 1.0.0: 12/nov/2010                                                   !
!  Coded by Saulo Freitas and Karla Longo                                       !
!  Contact: gmai@cptec.inpe.br - http://meioambiente.cptec.inpe.br              !
!###############################################################################!

MODULE util_geometry

IMPLICIT NONE

 TYPE pts_city

	REAL*4  :: y, x
	REAL*4  :: ang
	
 END TYPE pts_city

 TYPE T_mobileSrc
 
 	character(len=60)    		  :: cityName
  	character(len=60)  		  :: region
  	real		     		  :: area
  	real		     		  :: urbanArea
	real				  :: COEmission
	real				  :: NOXEmission
	real				  :: NMVOCEmission
  	real, dimension(2) 		  :: lowerBorder
  	real, dimension(2) 		  :: upperBorder
 	integer	     		          :: nUrbanVertexs
  	real, dimension(:), pointer   :: urbanXVertex
  	real, dimension(:), pointer   :: urbanYVertex

 END TYPE T_mobileSrc

TYPE(pts_city), ALLOCATABLE, DIMENSION(:) :: p_city
TYPE(pts_city)                            :: p_local
INTEGER                                   :: qte_p_city
INTEGER                                   :: max_x, max_y, min_x, min_y

CONTAINS

!****************************************************************************

SUBROUTINE PN_POLY(pol, qte, pt, ret)

 IMPLICIT NONE

 TYPE(pts_city), DIMENSION(:), INTENT(IN)  :: pol 
 TYPE(pts_city), INTENT(INOUT) :: pt
 INTEGER, INTENT(IN)         :: qte
 INTEGER, INTENT(OUT)        :: ret

 INTEGER :: i, cn
 REAL    :: vt

 cn = 0
 ret = 0
 
 DO i = 1, qte
 	IF( (pol(i)%y .LE. pt%y) .AND. (pol(i+1)%y .GT. pt%y) .OR. &
	    (pol(i)%y .GT. pt%y) .AND. (pol(i+1)%y .LE. pt%y))THEN
	    	vt = (pt%y - pol(i)%y)/(pol(i+1)%y - pol(i)%y)
		IF(pt%x .LT. (pol(i)%x + vt * (pol(i+1)%x - pol(i)%x))) cn = cn	+ 1
	ENDIF
 END DO
	   	
 ret = IAND(cn, 1)

 RETURN

END SUBROUTINE PN_POLY

!****************************************************************************
SUBROUTINE SET_ANCOR(pol, qte, max_x, max_y, min_x, min_y)

 IMPLICIT NONE
 
 TYPE(pts_city), DIMENSION(:), INTENT(INOUT) :: pol
 INTEGER, INTENT(INOUT)                      :: max_x, max_y, min_x, min_y
 INTEGER, INTENT(IN)           		     :: qte
 
 TYPE(pts_city) :: ancora
 INTEGER        :: i, j 
 
 ancora = pol(1)
 
 max_x = ancora%x
 min_x = ancora%x
 
 max_y = ancora%y
 min_y = ancora%y 
  
 DO i = 2, qte
 	IF(ancora%y .GT. pol(i)%y)THEN
		ancora = pol(i)
	ELSE
		IF(ancora%y .EQ. pol(i)%y)THEN
			IF(ancora%x .EQ. pol(i)%x) ancora = pol(i)
		END IF
	END IF
	
	IF(min_y .GT. pol(i)%y) min_y = pol(i)%y 
 	IF(max_y .LT. pol(i)%y) max_y = pol(i)%y
	IF(min_x .GT. pol(i)%x) min_x = pol(i)%x	
	IF(max_x .LT. pol(i)%x) max_x = pol(i)%x 
 END DO

 DO i = 1, qte
 	CALL CALC_THETA(ancora, pol(i))
 END DO

 pol(qte+1) = pol(1)

END SUBROUTINE SET_ANCOR

!****************************************************************************

SUBROUTINE CALC_THETA(P0, P1)

 IMPLICIT NONE

 TYPE(pts_city), INTENT(IN)  :: P0
 TYPE(pts_city), INTENT(INOUT) :: P1

 REAL*4 :: dx, dy, ax, ay, t

 dx = P1%x - P0%x
 dy = P1%y - P0%y

 ax = ABS(dx)
 ay = ABS(dy)

 IF(ax+ay == 0)THEN
         t = 0
 ELSE
         t = dy/(ax+ay)
 END IF

 IF(dx .LT. 0)THEN
         t = t -2
 ELSE
         t = t + 4
 END IF

 P1%ang = t * 90.

END SUBROUTINE CALC_THETA


SUBROUTINE loadMobileSrcMap(fileName, iMobile, voc_flag)
 
  !# parameter
  character(len=*),                         intent(in)  :: fileName
  type(T_mobileSrc), pointer, dimension(:) :: iMobile
  
  !# local
  integer :: nCities
  integer :: idx0
  integer :: idx1
  character :: dummyChar
  integer, optional :: voc_flag
  
  	open(222, file=fileName, status='old')
	read(222,*)nCities
	
	allocate(iMobile(nCities))
	
	do idx0 = 1, nCities
		read(222,*) iMobile(idx0)%cityName, dummyChar, iMobile(idx0)%region
		read(222,*) iMobile(idx0)%area, iMobile(idx0)%urbanArea
        if (present(voc_flag)) then
			read(222,*) iMobile(idx0)%COEmission, iMobile(idx0)%NOXEmission, iMobile(idx0)%NMVOCEmission
        else
			read(222,*) iMobile(idx0)%COEmission, iMobile(idx0)%NOXEmission
		end if
 		read(222,*) iMobile(idx0)%lowerBorder, iMobile(idx0)%upperBorder
		read(222,*) iMobile(idx0)%nUrbanVertexs
		
		print*, 'processing city: ', idx0, trim(iMobile(idx0)%cityName)
		
		allocate(iMobile(idx0)%urbanXVertex(iMobile(idx0)%nUrbanVertexs))
		allocate(iMobile(idx0)%urbanYVertex(iMobile(idx0)%nUrbanVertexs))
		
		do idx1 = 1, iMobile(idx0)%nUrbanVertexs
			read(222,*)iMobile(idx0)%urbanXVertex(idx1), iMobile(idx0)%urbanYVertex(idx1)
		end do
	
	end do
	
	close(222)
 
END SUBROUTINE loadMobileSrcMap

SUBROUTINE polyToModel(iMobileSrc,ng,n1,n2,xt,yt,deltax,deltay,plat,plon,rlon,rlat,grid_type, iCO, iNOX, iNMVOC)
 
 !# parameters
 type(T_mobileSrc), pointer, dimension(:)  :: iMobileSrc
 character(len=*),                         intent(in)  :: grid_type
 integer,                                  intent(in)  :: ng
 integer,                                  intent(in)  :: n1
 integer,                                  intent(in)  :: n2
 real, 					   intent(in)  :: deltax
 real, 					   intent(in)  :: deltay
 real, 					   intent(in)  :: plat
 real, 					   intent(in)  :: plon
 real, 					   intent(in)  :: xt(n1)
 real, 					   intent(in)  :: yt(n2)
 real, dimension(n1,n2), 		   intent(in)  :: rlon
 real, dimension(n1,n2), 		   intent(in)  :: rlat
 real, dimension(n1,n2),                   intent(out) :: iCO
 real, dimension(n1,n2),                   intent(out) :: iNOX
 real, dimension(n1,n2), optional,         intent(out) :: iNMVOC

 !# local
! real, parameter 		      :: resolX = 0.0083
! real, parameter 		      :: resolY = 0.0083
 real, parameter 		      :: resolX = 0.00021
 real, parameter 		      :: resolY = 0.00021
 integer 			      :: idx0
 integer 			      :: localX
 integer 			      :: localY
 integer 			      :: localNX
 integer 			      :: localNY
 integer 			      :: isInside
 integer 			      :: lat
 integer 			      :: lon
 integer 			      :: totalPixels
 real   			      :: localILat
 real    			      :: localILon
 real    			      :: sumCO
 real    			      :: sumNOX
 real    			      :: sumNMVOC 
 integer, allocatable, dimension(:,:) :: localGrid
 
 
 
 	allocate(localGrid(n1,n2))
	
	localGrid = 0.
 
 	do idx0 = 1, size(iMobileSrc)
	
		localILat = iMobileSrc(idx0)%upperBorder(2)
		localILon = iMobileSrc(idx0)%upperBorder(1)
		
		call update_emissions_by_city(ng,n1,n2,xt,yt,deltax,deltay,plat,plon,rlon,rlat, &
                          		      localILon,localILat,localX,localY,'none',grid_type)
		
		localILat = iMobileSrc(idx0)%lowerBorder(2)
		localILon = iMobileSrc(idx0)%lowerBorder(1)
		
		call update_emissions_by_city(ng,n1,n2,xt,yt,deltax,deltay,plat,plon,rlon,rlat, &
                          		      localILon,localILat,localX,localY,'none',grid_type)
		if(localX .eq. -999 .or. localY .eq. -999)then
		print *, 'pulei2', trim(iMobileSrc(idx0)%cityName)
		cycle
		end if
		
		localNX = abs(ceiling( (iMobileSrc(idx0)%upperBorder(1)-iMobileSrc(idx0)%lowerBorder(1))/resolX ))
		localNY = abs(ceiling( (iMobileSrc(idx0)%upperBorder(2)-iMobileSrc(idx0)%lowerBorder(2))/resolY ))
		
		totalPixels = 0
		localGrid   = 0.
		
		print*, 'city: ', idx0, trim(iMobileSrc(idx0)%cityName), ' - ', trim(iMobileSrc(idx0)%region)
		
		do lat = 1, localNY
			do lon = 1, localNX
				
				call insidePolygon(iMobileSrc(idx0)%urbanXVertex, iMobileSrc(idx0)%urbanYVertex, localILon, localILat, isInside)
				
				if(isInside .ne. 0)then
				
					call update_emissions_by_city(ng,n1,n2,xt,yt,deltax,deltay,plat,plon,rlon,rlat, &
                          					      localILon,localILat,localX,localY,'none',grid_type)
								
					if(localX .ne. -999 .and. localY .ne. -999)then
						
						localGrid(localX,localY) = 1 + localGrid(localX,localY)
				
						totalPixels = 1 + totalPixels
					end if
				
				end if
				
				localILon = localILon + resolX
			end do
			localILon = iMobileSrc(idx0)%lowerBorder(1)
			localILat = localILat + resolY
		end do
		
		print*, 'total pixels: ', totalPixels
		if(totalPixels .lt. 1) cycle
		
		sumCO  = 0.
		sumNOX = 0.
		sumNMVOC = 0.
		
		do lat = 1, n2
			do lon = 1, n1
				if(localGrid(lon,lat) .gt. 0.)then
					
					iCO(lon, lat)   = ( iMobileSrc(idx0)%COEmission * (localGrid(lon,lat)/(totalPixels*1.)) ) + iCO(lon, lat)
!					iCO(lon, lat)   = ( iMobileSrc(idx0)%COEmission * (localGrid(lon,lat)/(totalPixels -1.)) ) + iCO(lon, lat)
				sumCO  = iCO(lon, lat) + sumCO
					
					if(iCO(lon,lat) .lt. 0.)then
						print*, iMobileSrc(idx0)%COEmission, localGrid(lon,lat)
						stop  'emissao < 0'
					end if
					
					iNOX(lon, lat)  = ( (iMobileSrc(idx0)%NOXEmission)  * ((localGrid(lon,lat))/(totalPixels*1.)) ) + iNOX(lon, lat)
					sumNOX = iNOX(lon,lat) + sumNOX

                           if (present(iNMVOC)) then
					iNMVOC(lon, lat)  = ( (iMobileSrc(idx0)%NMVOCEmission)  * ((localGrid(lon,lat))/(totalPixels*1.)) ) + iNMVOC(lon, lat)
					sumNMVOC = iNMVOC(lon,lat) + sumNMVOC
                           end if
					
				end if
			
			
			end do
		end do
		
		
		! a soma das emissoes anteriores inviabiliza a conservacao de massa.
		
		!if( sumCO .ne. iMobileSrc(idx0)%COEmission) then
		!	print*, 'checking mass conservation>>> ',sumCO, iMobileSrc(idx0)%COEmission
		!else
		!	print*, 'checking mass conservation: ',sumCO, iMobileSrc(idx0)%COEmission
		!end if
		
		! Multiplicar emissoes anteriores por fator de correcao da atual?
		!if(sumCO  .ne. iMobileSrc(idx0)%COEmission)then
		!	iCO = iCO * (iMobileSrc(idx0)%COEmission/sumCO)
		!end if
					
		!if(sumNOX  .ne. iMobileSrc(idx0)%NOXEmission)then
		!	iNOX = iNOX * (iMobileSrc(idx0)%emissionNOX/sumNOX)
		!end if
					
	end do
	
	deallocate(localGrid)
  
 
 END SUBROUTINE polyToModel
 
 SUBROUTINE insidePolygon(xVertex, yVertex, pX, pY, isInside)
 
  	!# parameters
	real, dimension(:), intent(in)  :: xVertex
	real, dimension(:), intent(in)  :: yVertex
	real,               intent(in)  :: pX
	real,               intent(in)  :: pY
	integer,            intent(out) :: isInside
	
	!# local
	integer :: idx0
	real    :: dY
	
	
		isInside = 0
	
		if(size(xVertex) .ne. size(yVertex)) stop 'fatal error: subroutine insidePolygon - vertexs arrays are not conformable'
		
		do idx0 = 1, size(xVertex)-1
		
			if( (yVertex(idx0) .le. pY .and. yVertex(idx0+1) .gt. pY) .or. &
			    (yVertex(idx0) .gt. pY .and. yVertex(idx0+1) .le. pY)        )then
			    
				dY = (pY - yVertex(idx0)) / (yVertex(idx0+1) - yVertex(idx0)) 
				
				if( pX .lt. ( xVertex(idx0) + dY * (xVertex(idx0+1) - xVertex(idx0)) ) ) isInside = 1 + isInside 	    
			    
			end if
		
		end do
		
		isInside = iand(isInside,1)
		
 END SUBROUTINE insidePolygon

!****************************************************************************

END MODULE util_geometry
