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) :: minimum_distance
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)

cdr = pi / 180.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)) )
	allocate( Geographic_Coordinate_v3(dim_sizes(2),2) )
	status = sfrdata( sds_id, start, stride, edges, caliop_latitude )
		do i = 1 , dim_sizes(2)
		Geographic_Coordinate_v3(i,n) = caliop_latitude(1,i)  
		end do
	caliop_dim_sizes = dim_sizes(2)
	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 )
		do i = 1 , dim_sizes(2)
		Geographic_Coordinate_v3(i,n) = caliop_longitude(1,i)  
		end do
	end if

end do

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 )
		do i = 1 , dim_sizes(1)
			do j = 1 , dim_sizes(2)
			Geographic_Coordinate_mod03(i,j,k) = mod03_latitude(i,j)  
			end do
		end do
	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 )
		do i = 1 , dim_sizes(1)
			do j = 1 , dim_sizes(2)
			Geographic_Coordinate_mod03(i,j,k) = mod03_longitude(i,j)  
			end do
		end do
	end if


end do

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

Cell_Along_Swath_1km = dim_sizes(1) 
Cell_Across_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 )

!----------------------------------------------------------------------------------------!
! Data Screening (remove all modis pixels too far from Caliop Track)

nb_caliop_data_screen = 0

do i = 1, caliop_dim_sizes
	if( min_modis_latitude < caliop_latitude(1,i) .AND. &
	&   caliop_latitude(1,i) < max_modis_latitude )then
	if( min_modis_longitude < caliop_longitude(1,i) .AND. &
	&   caliop_longitude(1,i) < max_modis_longitude )then
		nb_caliop_data_screen = nb_caliop_data_screen + 1
	end if
	end if
end do

allocate( caliop_index_screen(nb_caliop_data_screen) )
allocate( caliop_latitude_screen(nb_caliop_data_screen) )
allocate( caliop_longitude_screen(nb_caliop_data_screen) )

nb_caliop_data_screen = 0

do i = 1, caliop_dim_sizes
	if( min_modis_latitude < caliop_latitude(1,i) .AND. &
	&   caliop_latitude(1,i) < max_modis_latitude )then
	if( min_modis_longitude < caliop_longitude(1,i) .AND. &
	&   caliop_longitude(1,i) < max_modis_longitude )then
		nb_caliop_data_screen = nb_caliop_data_screen + 1
		caliop_index_screen(nb_caliop_data_screen) = i 	
		caliop_latitude_screen(nb_caliop_data_screen) = caliop_latitude(1,i)
		caliop_longitude_screen(nb_caliop_data_screen) = caliop_longitude(1,i)
	end if
	end if
end do

min_caliop_longitude_screen = MINVAL( caliop_longitude_screen )
max_caliop_longitude_screen = MAXVAL( caliop_longitude_screen )
min_caliop_latitude_screen = MINVAL( caliop_latitude_screen )
max_caliop_latitude_screen = MAXVAL( caliop_latitude_screen )

!write(6,*) min_caliop_longitude_screen, max_caliop_longitude_screen
!write(6,*) min_caliop_latitude_screen, max_caliop_latitude_screen

nb_modis_data_screen = 0

do i = 1, Cell_Along_Swath_1km
	do j = 1, Cell_Across_Swath_1km
		if( min_caliop_latitude_screen < mod03_latitude(i,j) .AND. &
		&   mod03_latitude(i,j) < max_caliop_latitude_screen )then
		if( min_caliop_longitude_screen < mod03_longitude(i,j) .AND. & 
		&   mod03_longitude(i,j) < max_caliop_longitude_screen )then
		nb_modis_data_screen = nb_modis_data_screen + 1
		end if
		end if
	end do
end do

allocate( modis_index_screen(nb_modis_data_screen,2) )
allocate( modis_latitude_screen(nb_modis_data_screen) )
allocate( modis_longitude_screen(nb_modis_data_screen) )

nb_modis_data_screen = 0

do i = 1, Cell_Along_Swath_1km
	do j = 1, Cell_Across_Swath_1km
		if( min_caliop_latitude_screen < mod03_latitude(i,j) .AND. &
		&   mod03_latitude(i,j) < max_caliop_latitude_screen )then
		if( min_caliop_longitude_screen < mod03_longitude(i,j) .AND. &
		&   mod03_longitude(i,j) < max_caliop_longitude_screen )then
		nb_modis_data_screen = nb_modis_data_screen + 1
		modis_index_screen(nb_modis_data_screen,1) = i 
		modis_index_screen(nb_modis_data_screen,2) = j 
		modis_latitude_screen(nb_modis_data_screen) = mod03_latitude(i,j)
		modis_longitude_screen(nb_modis_data_screen) = mod03_longitude(i,j)
		end if
		end if
	end do
end do

min_modis_longitude_screen = MINVAL( modis_longitude_screen )
max_modis_longitude_screen = MAXVAL( modis_longitude_screen )
min_modis_latitude_screen = MINVAL( modis_latitude_screen )
max_modis_latitude_screen = MAXVAL( modis_latitude_screen )

!write(6,*) nb_modis_data_screen
!write(6,*) min_modis_longitude_screen, max_modis_longitude_screen
!write(6,*) min_modis_latitude_screen, max_modis_latitude_screen

!----------------------------------------------------------------------------------------!
! Modis Caliop Pixels Collocation 

!write(6,*) len( trim(myd03_file_name) )

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

do i = 1 , nb_caliop_data_screen
	minimum_distance = 1.0d0 ! Beyond, pixels are considered disjoints
	collocation = .false.
	do j = 1 , nb_modis_data_screen
		if( modis_latitude_screen(j) - 1.0 < caliop_latitude_screen(i) .AND. &
		&   caliop_latitude_screen(i)<modis_latitude_screen(j)+1)then
		if( modis_longitude_screen(j) - 1.0 < caliop_longitude_screen(i) .AND. &
		&   caliop_longitude_screen(i)<modis_longitude_screen(j)+1)then
			!---------- Ellipsoidal Earth projected to a plane ----------!
			diff_lat = caliop_latitude_screen(i) - modis_latitude_screen(j)
			diff_long = caliop_longitude_screen(i) - modis_longitude_screen(j)
			phim = ( modis_latitude_screen(j) + caliop_latitude_screen(i) ) / 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 < minimum_distance )then
				minimum_distance = distance
				min_j = j
				collocation = .true.
			end if
		end if
		end if
	end do
	if( collocation .eqv. .true. )then
	write(1,*) modis_index_screen(min_j,1), modis_index_screen(min_j,2), &
	&          caliop_index_screen(i)
	end if
end do

close(1)

end

About / Au sujet de

Ben
Created:
3 avril 2017 20:27:15

Updated:
3 avril 2017 20:27:15

License / Licence

MIT License

Abstract / Résumé