program colocated_modis_caliop

! Compilation (example)
!
! gfortran ModisCaliopCollocatedPixels.f90 
!          -I/hdf-4.2.4-snap4/hdf4/include 
!          -L/hdf-4.2.4-snap4/hdf4/lib -lmfhdf -ldf -ljpeg -lz 
!          -o ModisCaliopCollocatedPixels
!
! Execution (example)
!
! ModisCaliopCollocatedPixels CAL_LID_L2_01kmCLay_name MYD03_name
      
implicit none

integer :: i,j,k,l,n

integer, parameter :: pr = selected_real_kind(15,3)

integer :: ierr ! Read File Status: Success (if == 0) or Fail (if /= 0) 
integer :: uin=42 ! logical unita associated to the file 
character*300 caliop_file_name, myd03_file_name 
character*300 fn
character*2 day, heure, minute
integer :: offset

integer DFACC_READ , DFNT_INT32
parameter (DFACC_READ = 1, DFNT_INT32 = 24)
integer sfstart, sfselect, sfrdata, sfendacc, sfend 
integer sffattr , sfrattr, var_id, attr_id, err_code
real*4 :: scale_factors(20), offsets(20)
integer sd_id, sds_id, sds_index
integer status, n_attrs
integer sffinfo , sfginfo
integer n_datasets, n_file_attrs , index
integer rank, data_type
integer start(3), edges(3), stride(3)

integer X_LENGTH_01 , Y_LENGTH_01 
integer X_LENGTH_01_max , Y_LENGTH_01_max 
parameter ( X_LENGTH_01_max = 1354 , Y_LENGTH_01_max = 2040 ) 
integer :: nb_sds_index 
parameter ( nb_sds_index = 2 )
integer tab_sds_index(nb_sds_index)
integer MAX_NC_NAME , MAX_VAR_DIMS
parameter ( MAX_NC_NAME = 256 , MAX_VAR_DIMS = 32 )
integer dim_sizes( MAX_VAR_DIMS )
character*( MAX_NC_NAME ) name 

!----------------------------------------------------------------------------------------!
! Caliop Data 

real(kind=4), dimension(:,:), allocatable :: caliop_latitude ! data type 5 
real(kind=4), dimension(:,:), allocatable :: caliop_longitude ! data type 5 

real(kind=4), dimension(:,:), allocatable :: Geographic_Coordinate_v3 ! data type 5 

integer :: nb_sds_index_v3 
parameter ( nb_sds_index_v3 = 2 )
integer tab_sds_index_v3( nb_sds_index_v3 )      
data tab_sds_index_v3 /1,2/
integer :: caliop_dim_sizes

integer :: nb_caliop_data_screen
real(pr) :: min_caliop_longitude_screen, max_caliop_longitude_screen
real(pr) :: min_caliop_latitude_screen, max_caliop_latitude_screen
integer, dimension(:), allocatable :: caliop_index_screen
real(kind=4), dimension(:), allocatable :: caliop_latitude_screen ! data type 5
real(kind=4), dimension(:), allocatable :: caliop_longitude_screen ! data type 5 !

!----------------------------------------------------------------------------------------!
! MYD03 Data

real(kind=4), dimension(:,:), allocatable :: mod03_latitude ! data type 5
real(kind=4), dimension(:,:), allocatable :: mod03_longitude ! data type 5

real(kind=4), dimension(:,:,:), allocatable :: Geographic_Coordinate_mod03 ! data type 5

integer :: nb_sds_index_mod03
parameter ( nb_sds_index_mod03 = 2 )
integer tab_sds_index_mod03( nb_sds_index_mod03 )
data tab_sds_index_mod03 /0,1/
integer :: Cell_Along_Swath_1km , Cell_Across_Swath_1km

integer :: nb_modis_data_screen
real(pr) :: min_modis_longitude , max_modis_longitude
real(pr) :: min_modis_latitude , max_modis_latitude
real(pr) :: min_modis_longitude_screen, max_modis_longitude_screen
real(pr) :: min_modis_latitude_screen, max_modis_latitude_screen
integer, dimension(:,:), allocatable :: modis_index_screen
real(kind=4), dimension(:), allocatable :: modis_latitude_screen ! data type 5
real(kind=4), dimension(:), allocatable :: modis_longitude_screen ! data type 5

!----------------------------------------------------------------------------------------!
! Collocation Variables

real(pr) :: distance

real(pr) :: distance_pr



real(pr) :: minimum_distance


real(pr) :: distance_min

real(pr) :: earth_radius = 6371.0d0
real(pr) :: diff_lat, diff_long 
real(pr) :: phim
real(pr) :: K1, K2
real(pr) :: tmp_01, tmp_02, tmp_03
integer :: min_j, min_k

logical :: collocation

real(pr) :: cdr
real(pr) :: pi = dacos(-1.0d0)


integer :: i_pr, j_pr

integer :: i_min, j_min

logical :: find_min_distance
logical :: find_first_collocated_pixel
logical :: find_modis_pixel

integer :: tmp_n

cdr = pi / 180.0d0

minimum_distance = 1.0d0

!----------------------------------------------------------------------------------------!
! Read Caliop (Longitude & Latitude) 

call getarg(1,caliop_file_name)

sd_id = sfstart(caliop_file_name, DFACC_READ)
status = sffinfo( sd_id , n_datasets , n_file_attrs )

do n = 1 , nb_sds_index_v3

	sds_index = tab_sds_index_v3(n)
	sds_id = sfselect(sd_id, sds_index)
	status = sfginfo( sds_id , name , rank , dim_sizes , data_type , n_attrs )
	
	do i = 1 , rank
	start(i) = 0
	edges(i) = dim_sizes(i)
	stride(i) = 1
	end do

	!---------- caliop_latitude	----------!

	if( sds_index == 1 )then
	allocate( caliop_latitude(1,dim_sizes(2)) )
	status = sfrdata( sds_id, start, stride, edges, caliop_latitude )
	end if

	!---------- caliop_longitude ----------!

	if( sds_index == 2 )then
	allocate( caliop_longitude(1,dim_sizes(2)) )
	status = sfrdata( sds_id, start, stride, edges, caliop_longitude )
	end if

end do

caliop_dim_sizes = dim_sizes(2)

status = sfendacc(sds_id)
status = sfend(sd_id)

!----------------------------------------------------------------------------------------!
! Read MYD03 (Longitude & Latitude) 

call getarg(2,myd03_file_name)

sd_id = sfstart(myd03_file_name, DFACC_READ)
status = sffinfo( sd_id , n_datasets , n_file_attrs )

do n = 1 , nb_sds_index_mod03

	sds_index = tab_sds_index_mod03(n)
	sds_id = sfselect(sd_id, sds_index)
	status = sfginfo( sds_id , name , rank , dim_sizes , data_type , n_attrs )

	do i = 1 , rank
	start(i) = 0
	edges(i) = dim_sizes(i)
	stride(i) = 1
	end do

	!---------- myd03_longitude ----------!

	if( sds_index == 0 )then
	k = 1
	allocate( mod03_latitude(dim_sizes(1),dim_sizes(2)))
	allocate( Geographic_Coordinate_mod03(dim_sizes(1),dim_sizes(2),2) )
	status = sfrdata( sds_id, start, stride, edges, mod03_latitude )
	end if

	!---------- myd03_longitude ----------!
	
	if( sds_index == 1 )then
	k = 2
	allocate( mod03_longitude(dim_sizes(1),dim_sizes(2)))
	status = sfrdata( sds_id, start, stride, edges, mod03_longitude )
	end if


end do

status = sfendacc(sds_id)
status = sfend(sd_id)

Cell_Across_Swath_1km = dim_sizes(1)
Cell_Along_Swath_1km = dim_sizes(2) 

min_modis_longitude = MINVAL( mod03_longitude )
max_modis_longitude = MAXVAL( mod03_longitude )
min_modis_latitude = MINVAL( mod03_latitude )
max_modis_latitude = MAXVAL( mod03_latitude )

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

!write(6,*) 'min_modis_latitude, max_modis_latitude', min_modis_latitude, max_modis_latitude
!write(6,*) 'min_modis_longitude, max_modis_longitude', min_modis_longitude, max_modis_longitude

!write(6,*) shape(caliop_latitude), shape(mod03_longitude)
!write(6,*) Cell_Along_Swath_1km, Cell_Across_Swath_1km

!do i = 1, 10

!write(6,*) caliop_latitude(1,i), caliop_longitude(1,i)

!end do

!write(6,*)

!i = 650
!do j = 1, 10
!	write(6,*) mod03_latitude(i,j), mod03_longitude(i,j)
!end do
!write(6,*)

!----------------------------------------------------------------------------------------!
! Step 1: find the first modis-caliop collocated pixel

!open(1,file='calmod_index_new.txt')

offset = len( trim(myd03_file_name) ) - 41 ! myd03 file name has 41 letters

heure = myd03_file_name(offset+16:offset+17)
minute = myd03_file_name(offset+18:offset+19)

write(fn,'("calmod_index_    ")')

write(fn(14:15),'(a)') heure
write(fn(16:17),'(a)') minute

open(1,file=fn)

write(1,'(a)') caliop_file_name
write(1,'(a)') myd03_file_name

find_first_collocated_pixel = .False.

n = 1

do while( .not. find_first_collocated_pixel .and. n < caliop_dim_sizes )

	if( min_modis_latitude < caliop_latitude(1,n) .AND. &
	&   caliop_latitude(1,n) < max_modis_latitude )then
	if( min_modis_longitude < caliop_longitude(1,n) .AND. &
	&   caliop_longitude(1,n) < max_modis_longitude )then

		i_pr = 1
		j_pr = 1

		!--------------------------------------------------------------------------------!
		! Ellipsoidal Earth projected to a plane 
		
		diff_lat = caliop_latitude(1,n) - mod03_latitude(i_pr,j_pr)
		diff_long = caliop_longitude(1,n) - mod03_longitude(i_pr,j_pr)
	
		phim = ( mod03_latitude(i_pr,j_pr) + caliop_latitude(1,n) ) / 2.0
	
		K1 = 111.13209 - 0.56605 * cos( 2.0 * phim * cdr ) + &
		&    0.00120 * cos( 4.0 * phim * cdr )
	
		K2 = 111.41513 * cos( phim * cdr ) - 0.09455 * cos( 3.0 * phim * cdr ) + &
		&    0.00012 * cos( 5.0 * phim * cdr )
	
		tmp_01 = K1 * diff_lat
		tmp_02 = K2 * diff_long
	
		distance_pr = dsqrt( tmp_01**2 + tmp_02**2 )
		
		!--------------------------------------------------------------------------------!
	
		find_min_distance = .False.
		do while( .not. find_min_distance  )
				
		i = i_pr + 1
		j = j_pr
	
		!--------------------------------------------------------------------------------!
		! Ellipsoidal Earth projected to a plane 
		
		diff_lat = caliop_latitude(1,n) - mod03_latitude(i,j)
		diff_long = caliop_longitude(1,n) - mod03_longitude(i,j)
	
		phim = ( mod03_latitude(i,j) + caliop_latitude(1,n) ) / 2.0
	
		K1 = 111.13209 - 0.56605 * cos( 2.0 * phim * cdr ) + &
		&    0.00120 * cos( 4.0 * phim * cdr )
	
		K2 = 111.41513 * cos( phim * cdr ) - 0.09455 * cos( 3.0 * phim * cdr ) + &
		&    0.00012 * cos( 5.0 * phim * cdr )
	
		tmp_01 = K1 * diff_lat
		tmp_02 = K2 * diff_long
	
		distance = dsqrt( tmp_01**2 + tmp_02**2 )
		
		if( distance < distance_pr )then
			i_pr = i
			j_pr = j
			distance_pr = distance
			!write(6,*) i,j, distance	
		else
			find_min_distance = .True.
		end if
	
		!--------------------------------------------------------------------------------!
		
		end do

		!--------------------------------------------------------------------------------!	
	
		if( distance_pr < minimum_distance )then
			!minimum_distance = distance
			min_j = j
			collocation = .true.
			find_first_collocated_pixel = .True.				
		end if

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

	end if
	end if

n = n + 1

end do

!write(6,*) n-1,i_pr,j_pr, distance_pr 
!write(1,*) n-1,i_pr,j_pr, distance_pr 
write(1,*) i_pr,j_pr, n-1

!----------------------------------------------------------------------------------------!
! Step 2: iterate until no more collocated pixel can be found

tmp_n = 1

find_modis_pixel = .True.

do while( find_modis_pixel )

	!------------------------------------------------------------------------------------!
	! Ellipsoidal Earth projected to a plane 

	diff_lat = caliop_latitude(1,n) - mod03_latitude(i_pr,j_pr)
	diff_long = caliop_longitude(1,n) - mod03_longitude(i_pr,j_pr)

	phim = ( mod03_latitude(i_pr,j_pr) + caliop_latitude(1,n) ) / 2.0

	K1 = 111.13209 - 0.56605 * cos( 2.0 * phim * cdr ) + &
	&    0.00120 * cos( 4.0 * phim * cdr )

	K2 = 111.41513 * cos( phim * cdr ) - 0.09455 * cos( 3.0 * phim * cdr ) + &
	&    0.00012 * cos( 5.0 * phim * cdr )

	tmp_01 = K1 * diff_lat
	tmp_02 = K2 * diff_long

	distance_pr = dsqrt( tmp_01**2 + tmp_02**2 )

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

	find_min_distance = .False.

	do while( .not. find_min_distance  )

		distance_min = 999.9
		i_min = 1
		j_min = 1

		!--------------------------------------------------------------------------------!
		! neighbour 1

		i = i_pr + 1
		j = j_pr 

		if( i > 0 .and. i <=  Cell_Along_Swath_1km )then
		if( j > 0 .and. j <=  Cell_Along_Swath_1km )then

		diff_lat = caliop_latitude(1,n) - mod03_latitude(i,j)
		diff_long = caliop_longitude(1,n) - mod03_longitude(i,j)
	
		phim = ( mod03_latitude(i,j) + caliop_latitude(1,n) ) / 2.0
	
		K1 = 111.13209 - 0.56605 * cos( 2.0 * phim * cdr ) + &
		&    0.00120 * cos( 4.0 * phim * cdr )
	
		K2 = 111.41513 * cos( phim * cdr ) - 0.09455 * cos( 3.0 * phim * cdr ) + &
		&    0.00012 * cos( 5.0 * phim * cdr )
	
		tmp_01 = K1 * diff_lat
		tmp_02 = K2 * diff_long
	
		distance = dsqrt( tmp_01**2 + tmp_02**2 )

		if( distance < distance_min ) then
			distance_min = distance
			i_min = i
			j_min = j
		end if

		end if
		end if

		!--------------------------------------------------------------------------------!
		! neighbour 2

		i = i_pr + 1
		j = j_pr + 1

		if( i > 0 .and. i <=  Cell_Along_Swath_1km )then
		if( j > 0 .and. j <=  Cell_Along_Swath_1km )then
 
 		diff_lat = caliop_latitude(1,n) - mod03_latitude(i,j)
		diff_long = caliop_longitude(1,n) - mod03_longitude(i,j)
	
		phim = ( mod03_latitude(i,j) + caliop_latitude(1,n) ) / 2.0
	
		K1 = 111.13209 - 0.56605 * cos( 2.0 * phim * cdr ) + &
		&    0.00120 * cos( 4.0 * phim * cdr )
	
		K2 = 111.41513 * cos( phim * cdr ) - 0.09455 * cos( 3.0 * phim * cdr ) + &
		&    0.00012 * cos( 5.0 * phim * cdr )
	
		tmp_01 = K1 * diff_lat
		tmp_02 = K2 * diff_long
	
		distance = dsqrt( tmp_01**2 + tmp_02**2 )

		if( distance < distance_min ) then
			distance_min = distance
			i_min = i
			j_min = j
		end if
 
 		end if
		end if
 
		!--------------------------------------------------------------------------------!   
 		! neighbour 3
 		
		i = i_pr 
		j = j_pr + 1

		if( i > 0 .and. i <=  Cell_Along_Swath_1km )then
		if( j > 0 .and. j <=  Cell_Along_Swath_1km )then

		diff_lat = caliop_latitude(1,n) - mod03_latitude(i,j)
		diff_long = caliop_longitude(1,n) - mod03_longitude(i,j)
	
		phim = ( mod03_latitude(i,j) + caliop_latitude(1,n) ) / 2.0
	
		K1 = 111.13209 - 0.56605 * cos( 2.0 * phim * cdr ) + &
		&    0.00120 * cos( 4.0 * phim * cdr )
	
		K2 = 111.41513 * cos( phim * cdr ) - 0.09455 * cos( 3.0 * phim * cdr ) + &
		&    0.00012 * cos( 5.0 * phim * cdr )
	
		tmp_01 = K1 * diff_lat
		tmp_02 = K2 * diff_long
	
		distance = dsqrt( tmp_01**2 + tmp_02**2 )

		if( distance < distance_min ) then
			distance_min = distance
			i_min = i
			j_min = j
		end if

		end if
		end if

		!--------------------------------------------------------------------------------!
		! neighbour 4
		
		i = i_pr - 1
		j = j_pr + 1

		if( i > 0 .and. i <=  Cell_Along_Swath_1km )then
		if( j > 0 .and. j <=  Cell_Along_Swath_1km )then

		diff_lat = caliop_latitude(1,n) - mod03_latitude(i,j)
		diff_long = caliop_longitude(1,n) - mod03_longitude(i,j)
	
		phim = ( mod03_latitude(i,j) + caliop_latitude(1,n) ) / 2.0
	
		K1 = 111.13209 - 0.56605 * cos( 2.0 * phim * cdr ) + &
		&    0.00120 * cos( 4.0 * phim * cdr )
	
		K2 = 111.41513 * cos( phim * cdr ) - 0.09455 * cos( 3.0 * phim * cdr ) + &
		&    0.00012 * cos( 5.0 * phim * cdr )
	
		tmp_01 = K1 * diff_lat
		tmp_02 = K2 * diff_long
	
		distance = dsqrt( tmp_01**2 + tmp_02**2 )

		if( distance < distance_min ) then
			distance_min = distance
			i_min = i
			j_min = j
		end if

		end if
		end if

		!--------------------------------------------------------------------------------!
		! neighbour 5

		i = i_pr - 1 
		j = j_pr		

		if( i > 0 .and. i <=  Cell_Along_Swath_1km )then
		if( j > 0 .and. j <=  Cell_Along_Swath_1km )then

		diff_lat = caliop_latitude(1,n) - mod03_latitude(i,j)
		diff_long = caliop_longitude(1,n) - mod03_longitude(i,j)
	
		phim = ( mod03_latitude(i,j) + caliop_latitude(1,n) ) / 2.0
	
		K1 = 111.13209 - 0.56605 * cos( 2.0 * phim * cdr ) + &
		&    0.00120 * cos( 4.0 * phim * cdr )
	
		K2 = 111.41513 * cos( phim * cdr ) - 0.09455 * cos( 3.0 * phim * cdr ) + &
		&    0.00012 * cos( 5.0 * phim * cdr )
	
		tmp_01 = K1 * diff_lat
		tmp_02 = K2 * diff_long
	
		distance = dsqrt( tmp_01**2 + tmp_02**2 )

		if( distance < distance_min ) then
			distance_min = distance
			i_min = i
			j_min = j
		end if

		end if
		end if

		!--------------------------------------------------------------------------------!
		! neighbour 6

		i = i_pr - 1 
		j = j_pr - 1

		if( i > 0 .and. i <=  Cell_Along_Swath_1km )then
		if( j > 0 .and. j <=  Cell_Along_Swath_1km )then

		diff_lat = caliop_latitude(1,n) - mod03_latitude(i,j)
		diff_long = caliop_longitude(1,n) - mod03_longitude(i,j)
	
		phim = ( mod03_latitude(i,j) + caliop_latitude(1,n) ) / 2.0
	
		K1 = 111.13209 - 0.56605 * cos( 2.0 * phim * cdr ) + &
		&    0.00120 * cos( 4.0 * phim * cdr )
	
		K2 = 111.41513 * cos( phim * cdr ) - 0.09455 * cos( 3.0 * phim * cdr ) + &
		&    0.00012 * cos( 5.0 * phim * cdr )
	
		tmp_01 = K1 * diff_lat
		tmp_02 = K2 * diff_long
	
		distance = dsqrt( tmp_01**2 + tmp_02**2 )

		if( distance < distance_min ) then
			distance_min = distance
			i_min = i
			j_min = j
		end if

		end if
		end if

		!--------------------------------------------------------------------------------!
		! neighbour 7

		i = i_pr 
		j = j_pr - 1

		if( i > 0 .and. i <=  Cell_Along_Swath_1km )then
		if( j > 0 .and. j <=  Cell_Along_Swath_1km )then

		diff_lat = caliop_latitude(1,n) - mod03_latitude(i,j)
		diff_long = caliop_longitude(1,n) - mod03_longitude(i,j)
	
		phim = ( mod03_latitude(i,j) + caliop_latitude(1,n) ) / 2.0
	
		K1 = 111.13209 - 0.56605 * cos( 2.0 * phim * cdr ) + &
		&    0.00120 * cos( 4.0 * phim * cdr )
	
		K2 = 111.41513 * cos( phim * cdr ) - 0.09455 * cos( 3.0 * phim * cdr ) + &
		&    0.00012 * cos( 5.0 * phim * cdr )
	
		tmp_01 = K1 * diff_lat
		tmp_02 = K2 * diff_long
	
		distance = dsqrt( tmp_01**2 + tmp_02**2 )

		if( distance < distance_min ) then
			distance_min = distance
			i_min = i
			j_min = j
		end if

		end if
		end if

		!--------------------------------------------------------------------------------!
		! neighbour 8
		
		i = i_pr + 1
		j = j_pr - 1		

		if( i > 0 .and. i <=  Cell_Along_Swath_1km )then
		if( j > 0 .and. j <=  Cell_Along_Swath_1km )then

		diff_lat = caliop_latitude(1,n) - mod03_latitude(i,j)
		diff_long = caliop_longitude(1,n) - mod03_longitude(i,j)
	
		phim = ( mod03_latitude(i,j) + caliop_latitude(1,n) ) / 2.0
	
		K1 = 111.13209 - 0.56605 * cos( 2.0 * phim * cdr ) + &
		&    0.00120 * cos( 4.0 * phim * cdr )
	
		K2 = 111.41513 * cos( phim * cdr ) - 0.09455 * cos( 3.0 * phim * cdr ) + &
		&    0.00012 * cos( 5.0 * phim * cdr )
	
		tmp_01 = K1 * diff_lat
		tmp_02 = K2 * diff_long
	
		distance = dsqrt( tmp_01**2 + tmp_02**2 )

		if( distance < distance_min ) then
			distance_min = distance
			i_min = i
			j_min = j
		end if

		end if
		end if

		!write(6,*) distance_pr, i_pr, j_pr
		!write(6,*) distance_min, i_min, j_min
		!write(6,*) '-----------',  n, distance_min < distance_pr

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

		if( distance_min < distance_pr )then
			i_pr = i_min
			j_pr = j_min
			distance_pr = distance_min
		else
			find_min_distance = .True.
		end if

	end do

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

	if( distance_pr < minimum_distance )then
		!minimum_distance = distance
		min_j = j
		collocation = .true.
		find_first_collocated_pixel = .True.	
		!write(6,*) n,i_pr,j_pr, distance_pr, minimum_distance
		!write(1,*) n,i_pr,j_pr, distance_pr
		write(1,*) i_pr,j_pr,n
	end if

	if( distance_pr > 5.0 )then
		find_modis_pixel = .False.
	end if

	!tmp_n = tmp_n + 1
	!if( tmp_n > 10 ) STOP

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

n = n + 1

end do

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

close(1)

end

About / Au sujet de

Ben
Created:
3 avril 2017 20:59:24

Updated:
3 avril 2017 20:59:24

License / Licence

MIT License

Abstract / Résumé