Module SSMIS_Spatial_Average_Mod ! ! ! abstract: This routine reads BUFR format SSMIS radiance ! (brightness temperature) files, spatially ! averages the data using the AAPP averaging routines ! and writes the data back into BUFR format ! ! ! Program history log: ! 2011-11-18 collard - Original version (for ATMS) ! 2011-12-20 eliu - Modify to apply for SSMIS ! 2016-03-03 ejones - Add option for spatial averaging of GMI ! 2016-03-24 ejones - Add option for spatial averaging of AMSR2 ! use kinds, only: r_kind,r_double,i_kind use m_uniq use m_distance implicit none ! Declare module level parameters real(r_kind), parameter :: Missing_Value=1.e11_r_double CONTAINS SUBROUTINE SSMIS_Spatial_Average(BufrSat, Method, Num_Obs, NChanl, & FOV, Node_InOut, Time, Lat, Lon, BT_InOut, Error_status) IMPLICIT NONE ! Declare passed variables integer(i_kind) ,intent(in ) :: BufrSat integer(i_kind) ,intent(in ) :: Method ! 1=simple(1), 2=simple(2), 3=AAPP, 4=GMI (simple 1), 5=AMSR2 (simple1) integer(i_kind) ,intent(in ) :: Num_Obs, NChanl integer(i_kind) ,intent(in ) :: Fov(num_obs) integer(i_kind) ,intent(inout) :: Node_InOut(num_obs) real(r_kind) ,intent(in ) :: Time(Num_Obs) real(r_kind) ,intent(in ) :: Lat(Num_Obs) real(r_kind) ,intent(in ) :: Lon(Num_Obs) real(r_kind) ,intent(inout) :: BT_InOut(NChanl,Num_Obs) integer(i_kind) ,intent( out) :: Error_Status ! Declare local parameters ! integer(i_kind), parameter :: atms1c_h_wmosatid=224 integer(i_kind), parameter :: lninfile=15 ! integer(i_kind), parameter :: max_fov=96 integer(i_kind), parameter :: max_fov=60 integer(i_kind), parameter :: max_fov_gmi=221 integer(i_kind), parameter :: max_fov_amsr2=243 ! integer(i_kind), parameter :: max_obs=20000000 integer(i_kind), parameter :: as_node= 1_i_kind integer(i_kind), parameter :: ds_node=-1_i_kind real(r_kind), parameter :: btmin=70.0_r_kind real(r_kind), parameter :: btmax=320.0_r_kind real(r_kind), parameter :: sigma = 1.0_r_kind/25.0_r_kind ! real(r_kind), parameter :: sigma = 1.0_r_kind/50.0_r_kind ! Declare local variables character(100) :: infile character(30) :: Cline integer(i_kind) :: i,iscan,ifov,ichan,nchannels,wmosatid,version integer(i_kind) :: is,ip,ic integer(i_kind) :: iobs integer(i_kind) :: ios,max_scan,min_scan integer(i_kind) :: ns1,ns2,np1,np2 integer(i_kind) :: nscan integer(i_kind) :: ntime_scan integer(i_kind) :: delta_scan integer(i_kind) :: scanline_new integer(i_kind) :: nxaverage(nchanl),nyaverage(nchanl) integer(i_kind) :: channelnumber(nchanl) integer(i_kind) :: qc_distx(nchanl),qc_disty(nchanl) integer(i_kind), allocatable :: nodeinfo(:,:) integer(i_kind), allocatable :: scanline(:),scanline_back(:,:) real(r_kind), allocatable :: latitude(:,:), longitude(:,:) real(r_kind), allocatable :: time_scan(:) real(r_kind), allocatable :: dist_scan(:) real(r_kind), allocatable :: dist_scan_avg(:) real(r_kind), allocatable :: alat(:),alon(:),blat(:),blon(:) real(r_kind) :: dlat ! real(r_kind) :: time1,time2 real(r_kind) :: sampling_distx,sampling_disty,beamwidth(nchanl) real(r_kind) :: newwidth(nchanl),cutoff(nchanl) real(r_kind), allocatable, target :: bt_image(:,:,:) real(r_kind), allocatable :: bt_image_orig(:,:,:) real(r_kind), pointer :: bt_image1(:,:) real(r_kind) :: t1,t2,tdiff real(r_kind) :: lat1,lat2,lon1,lon2,dist,wgt real(r_kind) :: xnum,mta logical :: gaussian_wgt Error_Status=0 if (Method == 1) then ! simple averaging 1 gaussian_wgt = .false. ! write(*,*) 'SSMIS_Spatial_Average: using method from Banghua' write(*,*) 'SSMIS_Spatial_Average: bufrsat = ', BufrSat write(*,*) 'SSMIS_Spatial_Average: Gaussian Weighted Averaging = ', gaussian_wgt ! Determine scanline from time !============================== allocate(scanline(num_obs)) t1 = time(1) ! time for first scanline nscan = 1 ! first scanline scanline(1) = nscan do iobs = 2, num_obs t2 = time(iobs) tdiff = t2-t1 if (tdiff >= 0.00001_r_kind) then nscan = nscan+1 t1 = t2 endif scanline(iobs) = nscan enddo max_scan = maxval(scanline) write(*,*) 'SSMIS_Spatial_Average: max_scan,max_fov,nchanl = ', & max_scan,max_fov,nchanl ! Allocate and initialize variables allocate(bt_image_orig(max_fov,max_scan,nchanl)) allocate(latitude(max_fov,max_scan)) allocate(longitude(max_fov,max_scan)) allocate(nodeinfo(max_fov,max_scan)) allocate(scanline_back(max_fov,max_scan)) bt_image_orig(:,:,:)= 1000.0_r_kind latitude(:,:) = 1000.0_r_kind longitude(:,:) = 1000.0_r_kind scanline_back(:,:) = -1 nodeinfo(:,:) = 1000_i_kind ! Put data into 2D (fov vs. scanline) array do iobs = 1, num_obs latitude(fov(iobs),scanline(iobs)) = lat(iobs) longitude(fov(iobs),scanline(iobs)) = lon(iobs) bt_image_orig(fov(iobs),scanline(iobs),:)= bt_inout(:,iobs) scanline_back(fov(iobs),scanline(iobs)) = iobs enddo ! Determine AS/DS node information for each scanline loop1: do iscan = 1, max_scan-1 loop2: do ifov = 1, max_fov if (scanline_back(ifov,iscan) > 0 .and. scanline_back(ifov,iscan+1) > 0) then dlat = latitude(ifov,iscan+1)-latitude(ifov,iscan) if (dlat < 0.0_r_kind) then nodeinfo(:,iscan) = ds_node else nodeinfo(:,iscan) = as_node endif cycle loop1 endif enddo loop2 enddo loop1 nodeinfo(:,max_scan) = nodeinfo(:,max_scan-1) ! Do spatial averaging in the box centered on each fov for each channel !$omp parallel do schedule(dynamic,1)private(ic,iobs,iscan,ifov,ns1,ns2,np1,np2,xnum,mta,is,ip,lat1,lon1,lat2,lon2,dist,wgt) scan_loop: do iscan = 1, max_scan fov_loop: do ifov = 1, max_fov iobs = scanline_back(ifov,iscan) if (iobs >0) then node_inout(iobs) = nodeinfo(ifov,iscan) ! Define grid box (3 (scan direction) x 7 (satellite track dir)) ns1 = iscan-3 ns2 = iscan+3 if (ns1 < 1) ns1=1 if (ns2 > max_scan) ns2=max_scan np1 = ifov-1 np2 = ifov+1 if (np1 < 1) np1=1 if (np2 > max_fov) np2=max_fov channel_loop: do ic = 1, nchanl xnum = 0.0_r_kind mta = 0.0_r_kind if (any(bt_image_orig(np1:np2,ns1:ns2,ic) < btmin .or. & bt_image_orig(np1:np1,ns1:ns2,ic) > btmax)) then bt_inout(ic,iobs) = 1000.0_r_kind else ! Calculate distance of each fov to the center fov box_y1: do is = ns1, ns2 box_x1: do ip = np1, np2 lat1 = latitude(ifov,iscan) ! lat of the center fov lon1 = longitude(ifov,iscan) ! lon of the center fov lat2 = latitude(ip,is) lon2 = longitude(ip,is) dist = distance(lat1,lon1,lat2,lon2) if (dist > 100.0_r_kind) cycle box_x1 ! outside the box if (gaussian_wgt) then wgt = exp(-0.5_r_kind*(dist*sigma)*(dist*sigma)) else wgt = 1.0 endif xnum = xnum+wgt mta = mta +wgt*bt_image_orig(ip,is,ic) enddo box_x1 enddo box_y1 bt_inout(ic,iobs) = mta/xnum endif enddo channel_loop endif enddo fov_loop enddo scan_loop ! Deallocate arrays deallocate(nodeinfo,scanline,scanline_back) deallocate(latitude,longitude) deallocate(bt_image_orig) !============================================================================================================ ! Simple method 2 else if (Method == 2) then ! simple averaging 2 gaussian_wgt = .false. write(*,*) 'SSMIS_Spatial_Average: using method from Emily' write(*,*) 'SSMIS_Spatial_Average: bufrsat = ', BufrSat write(*,*) 'SSMIS_Spatial_Avearge: Gaussian Weighted Averaging = ', gaussian_wgt ! Determine scanline index from time ! Each uniq time is one uniq scanline ! call cpu_time(time1) ntime_scan = nuniq(time(1:num_obs)) allocate(scanline(num_obs)) allocate(time_scan(ntime_scan)) time_scan(1:ntime_scan) = time(uniq(time(1:num_obs),ntime_scan)) do iscan = 1, ntime_scan where((abs(time(1:num_obs)-time_scan(iscan)))<=0.0001_r_kind) scanline = iscan enddo max_scan = maxval(scanline) ! call cpu_time(time2) ! write(*,*)'CPU time for determining scanline index = ', time2-time1 ! write(*,*)'determine scanline index from time' ! write(*,*)'ntime_scan = ', ntime_scan ! write(*,*)'min/max scanline = ', minval(scanline), maxval(scanline) ! write(*,*)'min/max time = ', minval(time), maxval(time) ! write(*,*)'min/max time_scan = ', minval(time_scan), maxval(time_scan) ! Put lat/lon in 2D (fov,scanline) grid to calculate distance between scanlines ! call cpu_time(time1) allocate(latitude(max_fov,max_scan)) allocate(longitude(max_fov,max_scan)) allocate(scanline_back(max_fov,max_scan)) latitude(:,:) = 9999.0_r_kind ! filled value longitude(:,:) = 9999.0_r_kind ! filled value scanline_back(:,:) = -1 ! filled value do iobs = 1, num_obs latitude(fov(iobs),scanline(iobs)) = lat(iobs) longitude(fov(iobs),scanline(iobs)) = lon(iobs) scanline_back(fov(iobs),scanline(iobs))= iobs enddo ! Calculate distance between scanlines ! Reassign scanline according to distance between scanlines allocate(dist_scan(max_fov)) allocate(dist_scan_avg(max_scan)) dist_scan(:) = 0.0_r_kind dist_scan_avg(:) = 0.0_r_kind allocate(alat(max_fov),blat(max_fov),alon(max_fov),blon(max_fov)) scanline_new = 1 do iscan = 2, max_scan ! loop through scanlines dist_scan(:) = 0.0_r_kind alat(:) = latitude(:,iscan-1) alon(:) = longitude(:,iscan-1) blat(:) = latitude(:,iscan) blon(:) = longitude(:,iscan) ! where(alat>=1000.0_r_kind .or. alon>=1000.0_r_kind .or. blat>=1000.0_r_kind .or. blon>=1000.0_r_kind) ! missing FOVs where(scanline_back(:,iscan) <= 0 .or. scanline_back(:,iscan-1) <= 0) ! missing FOVs alat(:)=0.0_r_kind blat(:)=0.0_r_kind alon(:)=0.0_r_kind blon(:)=0.0_r_kind end where dist_scan = distance(alat,alon,blat,blon) dist_scan_avg(iscan) = sum(dist_scan)/count(dist_scan>0.0000001_r_kind) delta_scan = nint(dist_scan_avg(iscan)/12.5_r_kind) scanline_new = scanline_new+delta_scan do ifov = 1, max_fov iobs = scanline_back(ifov,iscan) if (iobs > 0) scanline(iobs) = scanline_new enddo enddo max_scan = maxval(scanline) min_scan = minval(scanline) ! call cpu_time(time2) ! write(*,*) 'CPU time for determining dist and reassign scanlines = ', time2-time1 ! write(*,*) 'SSMIS_Spatial_Average: reassigned max_scan = ', max_scan ! write(*,*) 'SSMIS_Spatial_Average: reassigned max_fov = ', max_fov ! write(*,*) 'SSMIS_Spatial_Average: reassigned nchanl = ', nchanl deallocate(time_scan,dist_scan,dist_scan_avg) deallocate(latitude,longitude,scanline_back) deallocate(alat,alon,blat,blon) ! Allocate and initialize variables allocate(latitude(max_fov,max_scan)) allocate(longitude(max_fov,max_scan)) allocate(bt_image(max_fov,max_scan,nchanl)) allocate(bt_image_orig(max_fov,max_scan,nchanl)) allocate(nodeinfo(max_fov,max_scan)) allocate(scanline_back(max_fov,max_scan)) bt_image(:,:,:) = 1000.0_r_kind bt_image_orig(:,:,:)= 1000.0_r_kind latitude(:,:) = 1000.0_r_kind longitude(:,:) = 1000.0_r_kind scanline_back(:,:) = -1 nodeinfo(:,:) = 1000_i_kind ! Put data in 2D (fov,scanline) grid to perform noise reduction do i=1,num_obs latitude(fov(i),scanline(i)) = lat(i) longitude(fov(i),scanline(i)) = lon(i) bt_image(fov(i),scanline(i),:) = bt_inout(:,i) bt_image_orig(fov(i),scanline(i),:)= bt_inout(:,i) scanline_back(fov(i),scanline(i)) = i enddo ! Determine AS/DS node information for each scanline loop3: do iscan = 1, max_scan-1 loop4: do ifov = 1, max_fov if (scanline_back(ifov,iscan) > 0 .and. scanline_back(ifov,iscan+1) > 0) then dlat = latitude(ifov,iscan+1)-latitude(ifov,iscan) if (dlat < 0.0_r_kind) then nodeinfo(:,iscan) = ds_node else nodeinfo(:,iscan) = as_node endif cycle loop3 endif enddo loop4 enddo loop3 nodeinfo(:,max_scan) = nodeinfo(:,max_scan-1) ! write(*,*) 'ssmis_spatial_average' ! write(*,*) 'min/max fov = ', minval(fov), maxval(fov) ! write(*,*) 'min/max time = ', minval(time), maxval(time) ! write(*,*) 'min/max scanline_back = ', minval(scanline_back), maxval(scanline_back) ! write(*,*) 'min/max scanline = ', minval(scanline), maxval(scanline) ! write(*,*) 'min/max scan = ', minval(scan), maxval(scan) ! Do spatial averaging in the box centered on each fov for each channel write(*,*) 'SSMIS_Spatial_Average: do spatial averaging for noise reduction ' scan_loop2: do iscan = 1, max_scan fov_loop2: do ifov = 1, max_fov ! Define grid box (3 (scan direction) x 7 (satellite track dir)) ns1 = iscan-3 ns2 = iscan+3 if (ns1 < 1) ns1=1 if (ns2 > max_scan) ns2=max_scan np1 = ifov-1 np2 = ifov+1 if (np1 < 1) np1=1 if (np2 > max_fov) np2=max_fov channel_loop2: do ic = 1, nchanl xnum = 0.0_r_kind mta = 0.0_r_kind if (any(bt_image_orig(np1:np2,ns1:ns2,ic) < btmin .or. bt_image_orig(np1:np1,ns1:ns2,ic) > btmax)) then bt_image(ifov,iscan,ic) = 1000.0_r_kind else ! Calculate distance of each fov to the center fov box_y2: do is = ns1, ns2 box_x2: do ip = np1, np2 lat1 = latitude(ifov,iscan) ! lat of the center fov lon1 = longitude(ifov,iscan) ! lon of the center fov lat2 = latitude(ip,is) lon2 = longitude(ip,is) dist = distance(lat1,lon1,lat2,lon2) if (dist > 100.0_r_kind) cycle box_x2 ! outside the box if (gaussian_wgt) then wgt = exp(-0.5_r_kind*(dist*sigma)*(dist*sigma)) else wgt = 1.0 endif xnum = xnum+wgt mta = mta +wgt*bt_image_orig(ip,is,ic) enddo box_x2 enddo box_y2 bt_image(ifov,iscan,ic) = mta/xnum endif enddo channel_loop2 enddo fov_loop2 enddo scan_loop2 do iscan = 1, max_scan do ifov = 1, max_fov if (scanline_back(ifov,iscan) >0) then bt_inout(:,scanline_back(ifov,iscan)) = bt_image(ifov,iscan,:) node_inout(scanline_back(ifov,iscan)) = nodeinfo(ifov,iscan) endif enddo enddo deallocate(nodeinfo,scanline,scanline_back) deallocate(latitude,longitude) deallocate(bt_image_orig,bt_image) !============================================================================================================ else if (Method == 3) then ! AAPP method write(*,*) 'SSMIS_Spatial_Average: using AAPP method' if (bufrsat == 249) infile= 'ssmis_f16_beamwidth.txt' if (bufrsat == 285) infile= 'ssmis_f17_beamwidth.txt' if (bufrsat == 286) infile= 'ssmis_f18_beamwidth.txt' if (bufrsat == 287) infile= 'ssmis_f19_beamwidth.txt' ! Read the beamwidth requirements OPEN(lninfile,file=infile,form='formatted',status='old', & iostat=ios) IF (ios /= 0) THEN WRITE(*,*) 'Unable to open ',trim(infile) Error_Status=1 RETURN ENDIF wmosatid=999 read(lninfile,'(a30)',iostat=ios) Cline DO WHILE (wmosatid .NE. BufrSat .AND. ios == 0) DO WHILE (Cline(1:1) == '#') read(lninfile,'(a30)') Cline ENDDO READ(Cline,*) wmosatid read(lninfile,'(a30)') Cline DO WHILE (Cline(1:1) == '#') read(lninfile,'(a30)') Cline ENDDO READ(Cline,*) version read(lninfile,'(a30)') Cline DO WHILE (Cline(1:1) == '#') read(lninfile,'(a30)') Cline ENDDO READ(Cline,*) sampling_distx, sampling_disty read(lninfile,'(a30)') Cline DO WHILE (Cline(1:1) == '#') read(lninfile,'(a30)') Cline ENDDO READ(Cline,*) nchannels read(lninfile,'(a30)') Cline if (nchannels > 0) then DO ichan=1,nchannels read(lninfile,'(a30)') Cline DO WHILE (Cline(1:1) == '#') read(lninfile,'(a30)') Cline ENDDO READ(Cline,*) channelnumber(ichan),beamwidth(ichan), & & newwidth(ichan),cutoff(ichan),nxaverage(ichan), & & nyaverage(ichan), qc_distx(ichan), qc_disty(ichan) ENDDO end if read(lninfile,'(a30)',iostat=ios) Cline ENDDO IF (wmosatid /= BufrSat) THEN WRITE(*,*) 'SSMIS_Spatial_Averaging: sat id not matched in ', trim(infile) Error_Status=1 RETURN ENDIF CLOSE(lninfile) !orig for ATMS ! ! Determine scanline from time ! MinTime = MINVAL(Time) ! ALLOCATE(Scanline(Num_Obs)) ! Scanline(:) = NINT((Time(1:Num_Obs)-MinTime)/Scan_Interval)+1 ! Max_Scan=MAXVAL(Scanline) ! Min_Scan=MINVAL(Scanline) ! Determine scanline index from time ! Each uniq time is one uniq scanline ! call cpu_time(time1) ntime_scan = nuniq(time(1:num_obs)) allocate(scanline(num_obs)) allocate(time_scan(ntime_scan)) time_scan(1:ntime_scan) = time(uniq(time(1:num_obs),ntime_scan)) do iscan = 1, ntime_scan where((abs(time(1:num_obs)-time_scan(iscan)))<=0.0001_r_kind) scanline = iscan enddo max_scan = maxval(scanline) ! call cpu_time(time2) ! write(*,*)'CPU time for determining scanline index = ', time2-time1 ! write(*,*)'determine scanline index from time' ! write(*,*)'ntime_scan = ', ntime_scan ! write(*,*)'min/max scanline = ', minval(scanline), maxval(scanline) ! write(*,*)'min/max time = ', minval(time), maxval(time) ! write(*,*)'min/max time_scan = ', minval(time_scan), maxval(time_scan) ! Put lat/lon in 2D (fov,scanline) grid to calculate distance between scanlines ! call cpu_time(time1) allocate(latitude(max_fov,max_scan)) allocate(longitude(max_fov,max_scan)) allocate(scanline_back(max_fov,max_scan)) latitude(:,:) = 9999.0_r_kind ! filled value longitude(:,:) = 9999.0_r_kind ! filled value scanline_back(:,:) = -1 ! filled value do iobs = 1, num_obs latitude(fov(iobs),scanline(iobs)) = lat(iobs) longitude(fov(iobs),scanline(iobs)) = lon(iobs) scanline_back(fov(iobs),scanline(iobs))= iobs enddo ! Calculate distance between scanlines ! Reassign scanline according to distance between scanlines allocate(dist_scan(max_fov)) allocate(dist_scan_avg(max_scan)) dist_scan(:) = 0.0_r_kind dist_scan_avg(:) = 0.0_r_kind allocate(alat(max_fov),blat(max_fov),alon(max_fov),blon(max_fov)) scanline_new = 1 do iscan = 2, max_scan ! loop through scanlines dist_scan(:) = 0.0_r_kind alat(:) = latitude(:,iscan-1) alon(:) = longitude(:,iscan-1) blat(:) = latitude(:,iscan) blon(:) = longitude(:,iscan) ! where(alat>=1000.0_r_kind .or. alon>=1000.0_r_kind .or. blat>=1000.0_r_kind .or. blon>=1000.0_r_kind) ! missing FOVs where(scanline_back(:,iscan) <= 0 .or. scanline_back(:,iscan-1) <= 0) ! missing FOVs alat(:)=0.0_r_kind blat(:)=0.0_r_kind alon(:)=0.0_r_kind blon(:)=0.0_r_kind end where dist_scan = distance(alat,alon,blat,blon) dist_scan_avg(iscan) = sum(dist_scan)/count(dist_scan>0.0000001_r_kind) delta_scan = nint(dist_scan_avg(iscan)/12.5_r_kind) scanline_new = scanline_new+delta_scan do ifov = 1, max_fov iobs = scanline_back(ifov,iscan) if (iobs > 0) scanline(iobs) = scanline_new enddo enddo max_scan = maxval(scanline) min_scan = minval(scanline) ! call cpu_time(time2) ! write(*,*)'CPU time for determining dist and reassign scanlines = ', time2-time1 ! write(*,*) 'SSMIS_Spatial_Average: reassigned max_scan = ', max_scan ! write(*,*) 'SSMIS_Spatial_Average: reassigned max_fov = ', max_fov ! write(*,*) 'SSMIS_Spatial_Average: reassigned nchanl = ', nchanl deallocate(time_scan,dist_scan,dist_scan_avg) deallocate(latitude,longitude,scanline_back) ! Put lat/lon in 2D (fov,scanline) grid to perform noise reduction ! call cpu_time(time1) ALLOCATE(BT_Image(Max_FOV,Max_Scan,nchanl)) ALLOCATE(BT_Image_Orig(Max_FOV,Max_Scan,nchanl)) ALLOCATE(Scanline_Back(Max_FOV,Max_Scan)) BT_Image(:,:,:) = 1000.0_r_kind ScanLine_Back(:,:) = -1 DO I=1,Num_Obs bt_image(FOV(I),Scanline(I),:)=bt_inout(:,I) bt_image_orig(FOV(I),Scanline(I),:)=bt_inout(:,I) Scanline_Back(FOV(I),Scanline(I))=I END DO ! write(*,*) 'ssmis_spatial_average' ! write(*,*) 'min/max fov = ', minval(fov), maxval(fov) ! write(*,*) 'min/max time = ', minval(time), maxval(time) ! write(*,*) 'min/max scanline_back = ', minval(scanline_back), maxval(scanline_back) ! write(*,*) 'min/max scanline = ', minval(scanline), maxval(scanline) ! write(*,*) 'min/max scan = ', minval(scan), maxval(scan) ! write(*,*) 'min/max bt_inout ', minval(bt_inout), maxval(bt_inout) ! Noise reduction using FFT DO IChan=1,nchanl bt_image1 => bt_image(:,:,ichan) ! If the channel number is present in the channelnumber array we should process it ! (otherwise bt_inout just keeps the same value): IF (ANY(channelnumber(1:nchannels) == ichan)) THEN CALL MODIFY_BEAMWIDTH ( max_fov, max_scan, bt_image1, & sampling_distx, sampling_disty, beamwidth(ichan), newwidth(ichan), & cutoff(ichan), nxaverage(ichan), nyaverage(ichan), & qc_distx(ichan), qc_disty(ichan), IOS) IF (IOS == 0) THEN do iscan=1,max_scan do ifov=1,max_fov if (Scanline_Back(IFov, IScan) > 0) & bt_inout(channelnumber(ichan),Scanline_Back(IFov, IScan)) = & BT_Image1(ifov,iscan) end do end do ELSE Error_Status=1 RETURN END IF END IF END DO DEALLOCATE(BT_Image, Scanline, Scanline_Back) NULLIFY(BT_Image1) ! call cpu_time(time2) ! write(*,*)'CPU time for noise reduction = ', time2-time1 endif ! method=3 !============================================================================================================ ! Simple method for GMI (like method 1) if (Method == 4) then ! simple averaging 1 gaussian_wgt = .false. ! write(*,*) 'SSMIS_Spatial_Average for GMI: using method from Banghua' write(*,*) 'SSMIS_Spatial_Average for GMI: bufrsat = ', BufrSat write(*,*) 'SSMIS_Spatial_Average for GMI: Gaussian Weighted Averaging =',gaussian_wgt ! Determine scanline from time !============================== allocate(scanline(num_obs)) t1 = time(1) ! time for first scanline nscan = 1 ! first scanline scanline(1) = nscan do iobs = 2, num_obs t2 = time(iobs) tdiff = t2-t1 if (tdiff >= 0.00001_r_kind) then nscan = nscan+1 t1 = t2 endif scanline(iobs) = nscan enddo max_scan = maxval(scanline) write(*,*) 'SSMIS_Spatial_Average for GMI:max_scan,max_fov,nchanl = ', & max_scan,max_fov,nchanl ! Allocate and initialize variables allocate(bt_image_orig(max_fov_gmi,max_scan,nchanl)) allocate(latitude(max_fov_gmi,max_scan)) allocate(longitude(max_fov_gmi,max_scan)) allocate(nodeinfo(max_fov_gmi,max_scan)) allocate(scanline_back(max_fov_gmi,max_scan)) bt_image_orig(:,:,:)= 1000.0_r_kind latitude(:,:) = 1000.0_r_kind longitude(:,:) = 1000.0_r_kind scanline_back(:,:) = -1 nodeinfo(:,:) = 1000_i_kind ! Put data into 2D (fov vs. scanline) array do iobs = 1, num_obs latitude(fov(iobs),scanline(iobs)) = lat(iobs) longitude(fov(iobs),scanline(iobs)) = lon(iobs) bt_image_orig(fov(iobs),scanline(iobs),:)= bt_inout(:,iobs) scanline_back(fov(iobs),scanline(iobs)) = iobs enddo ! Determine AS/DS node information for each scanline gmi_loop1: do iscan = 1, max_scan-1 gmi_loop2: do ifov = 1, max_fov_gmi if (scanline_back(ifov,iscan) > 0 .and. scanline_back(ifov,iscan+1)> 0) then dlat = latitude(ifov,iscan+1)-latitude(ifov,iscan) if (dlat < 0.0_r_kind) then nodeinfo(:,iscan) = ds_node else nodeinfo(:,iscan) = as_node endif cycle gmi_loop1 endif enddo gmi_loop2 enddo gmi_loop1 nodeinfo(:,max_scan) = nodeinfo(:,max_scan-1) ! Do spatial averaging in the box centered on each fov for each channel gmi_scan_loop: do iscan = 1, max_scan gmi_fov_loop: do ifov = 1, max_fov_gmi iobs = scanline_back(ifov,iscan) if (iobs >0) then node_inout(iobs) = nodeinfo(ifov,iscan) gmi_channel_loop: do ic = 1, nchanl ! Define grid box by channel - ! Ch 1-2: 1 scan direction, 1 track direction ! Ch 3-13: 3 scan direction, 3 track direction if ((ic == 1) .or. (ic == 2)) then ns1 = iscan ns2 = iscan if (ns1 < 1) ns1=1 if (ns2 > max_scan) ns2=max_scan np1 = ifov np2 = ifov if (np1 < 1) np1=1 if (np2 > max_fov_gmi) np2=max_fov_gmi else if ((ic > 2) .and. (ic < 14)) then ns1 = iscan-1 ns2 = iscan+1 if (ns1 < 1) ns1=1 if (ns2 > max_scan) ns2=max_scan np1 = ifov-1 np2 = ifov+1 if (np1 < 1) np1=1 if (np2 > max_fov_gmi) np2=max_fov_gmi endif xnum = 0.0_r_kind mta = 0.0_r_kind if (any(bt_image_orig(np1:np2,ns1:ns2,ic) < btmin .or. & bt_image_orig(np1:np1,ns1:ns2,ic) > btmax)) then bt_inout(ic,iobs) = 1000.0_r_kind else ! Calculate distance of each fov to the center fov gmi_box_y1: do is = ns1, ns2 gmi_box_x1: do ip = np1, np2 lat1 = latitude(ifov,iscan) ! lat of the center fov lon1 = longitude(ifov,iscan) ! lon of the center fov lat2 = latitude(ip,is) lon2 = longitude(ip,is) dist = distance(lat1,lon1,lat2,lon2) if (dist > 50.0_r_kind) cycle gmi_box_x1 ! outside the box if (gaussian_wgt) then wgt = exp(-0.5_r_kind*(dist/sigma)*(dist/sigma)) else wgt = 1.0 endif xnum = xnum+wgt mta = mta +wgt*bt_image_orig(ip,is,ic) enddo gmi_box_x1 enddo gmi_box_y1 bt_inout(ic,iobs) = mta/xnum endif enddo gmi_channel_loop endif enddo gmi_fov_loop enddo gmi_scan_loop ! Deallocate arrays deallocate(nodeinfo,scanline,scanline_back) deallocate(latitude,longitude) deallocate(bt_image_orig) endif ! Method=4 !============================================================================================================ ! Simple method for AMSR2 (like method 1) if (Method == 5) then ! simple averaging 1 gaussian_wgt = .false. ! write(*,*) 'SSMIS_Spatial_Average for AMSR2: using method from Banghua' write(*,*) 'SSMIS_Spatial_Average for AMSR2: bufrsat = ', BufrSat write(*,*) 'SSMIS_Spatial_Average for AMSR2: Gaussian Weighted Averaging=',gaussian_wgt ! Determine scanline from time !============================== allocate(scanline(num_obs)) t1 = time(1) ! time for first scanline nscan = 1 ! first scanline scanline(1) = nscan do iobs = 2, num_obs t2 = time(iobs) tdiff = t2-t1 if (tdiff >= 0.00001_r_kind) then nscan = nscan+1 t1 = t2 endif scanline(iobs) = nscan enddo max_scan = maxval(scanline) write(*,*) 'SSMIS_Spatial_Average for AMSR2:max_scan,max_fov,nchanl = ', & max_scan,max_fov,nchanl ! Allocate and initialize variables allocate(bt_image_orig(max_fov_amsr2,max_scan,nchanl)) allocate(latitude(max_fov_amsr2,max_scan)) allocate(longitude(max_fov_amsr2,max_scan)) allocate(nodeinfo(max_fov_amsr2,max_scan)) allocate(scanline_back(max_fov_amsr2,max_scan)) bt_image_orig(:,:,:)= 1000.0_r_kind latitude(:,:) = 1000.0_r_kind longitude(:,:) = 1000.0_r_kind scanline_back(:,:) = -1 nodeinfo(:,:) = 1000_i_kind ! Put data into 2D (fov vs. scanline) array do iobs = 1, num_obs latitude(fov(iobs),scanline(iobs)) = lat(iobs) longitude(fov(iobs),scanline(iobs)) = lon(iobs) bt_image_orig(fov(iobs),scanline(iobs),:)= bt_inout(:,iobs) scanline_back(fov(iobs),scanline(iobs)) = iobs enddo ! Determine AS/DS node information for each scanline amsr2_loop1: do iscan = 1, max_scan-1 amsr2_loop2: do ifov = 1, max_fov_amsr2 if (scanline_back(ifov,iscan) > 0 .and. scanline_back(ifov,iscan+1)> 0) then dlat = latitude(ifov,iscan+1)-latitude(ifov,iscan) if (dlat < 0.0_r_kind) then nodeinfo(:,iscan) = ds_node else nodeinfo(:,iscan) = as_node endif cycle amsr2_loop1 endif enddo amsr2_loop2 enddo amsr2_loop1 nodeinfo(:,max_scan) = nodeinfo(:,max_scan-1) ! Do spatial averaging in the box centered on each fov for each channel amsr2_scan_loop: do iscan = 1, max_scan amsr2_fov_loop: do ifov = 1, max_fov_amsr2 iobs = scanline_back(ifov,iscan) if (iobs >0) then node_inout(iobs) = nodeinfo(ifov,iscan) amsr2_channel_loop: do ic = 1, nchanl ! Define grid box by channel - ! Ch 1-6: 1 scan direction, 1 track direction ! Ch 7-14: 3 scan direction, 3 track direction if ((ic >= 1) .and. (ic <= 6)) then ns1 = iscan ns2 = iscan if (ns1 < 1) ns1=1 if (ns2 > max_scan) ns2=max_scan np1 = ifov np2 = ifov if (np1 < 1) np1=1 if (np2 > max_fov_gmi) np2=max_fov_amsr2 else if ((ic >= 7) .and. (ic <= 14)) then ns1 = iscan-1 ns2 = iscan+1 if (ns1 < 1) ns1=1 if (ns2 > max_scan) ns2=max_scan np1 = ifov-1 np2 = ifov+1 if (np1 < 1) np1=1 if (np2 > max_fov_amsr2) np2=max_fov_amsr2 endif xnum = 0.0_r_kind mta = 0.0_r_kind if (any(bt_image_orig(np1:np2,ns1:ns2,ic) < btmin .or. & bt_image_orig(np1:np1,ns1:ns2,ic) > btmax)) then bt_inout(ic,iobs) = 1000.0_r_kind else ! Calculate distance of each fov to the center fov amsr2_box_y1: do is = ns1, ns2 amsr2_box_x1: do ip = np1, np2 lat1 = latitude(ifov,iscan) ! lat of the center fov lon1 = longitude(ifov,iscan) ! lon of the center fov lat2 = latitude(ip,is) lon2 = longitude(ip,is) dist = distance(lat1,lon1,lat2,lon2) if (dist > 50.0_r_kind) cycle amsr2_box_x1 ! outside the box if (gaussian_wgt) then wgt = exp(-0.5_r_kind*(dist/sigma)*(dist/sigma)) else wgt = 1.0 endif xnum = xnum+wgt mta = mta +wgt*bt_image_orig(ip,is,ic) enddo amsr2_box_x1 enddo amsr2_box_y1 bt_inout(ic,iobs) = mta/xnum endif enddo amsr2_channel_loop endif enddo amsr2_fov_loop enddo amsr2_scan_loop ! Deallocate arrays deallocate(nodeinfo,scanline,scanline_back) deallocate(latitude,longitude) deallocate(bt_image_orig) endif ! Method=5 END Subroutine SSMIS_Spatial_Average SUBROUTINE MODIFY_BEAMWIDTH ( nx, ny, image, sampling_distx, sampling_disty, & beamwidth, newwidth, mtfcutoff, nxaverage, nyaverage, qc_distx, qc_disty, & Error) !----------------------------------------- ! Name: $Id$ ! ! Purpose: ! Manipulate the effective beam width of an image. For example, convert ATMS ! to AMSU-A-like resolution while reducing the noise. ! ! Method: ! 1) Pad the image to a power of 2 in each dimension. ! If FFT technique is to be used then: ! 2) Assuming Gaussian beam shapes, calcluate the input and output Modulation ! Transfer Functions (MTF). ! 3) FFT image to frequency domain (2-D). ! 4) Multiply by output MTF divided by input MTF. If a cut-off is specified ! (when attempting to make the beam width narrower), attenuate further ! by an exponential function - factor of 2 at the cutoff. ! 5) FFT back to image domain ! Finally, ! 6) Over-write the input image, with averaging if requested. ! ! COPYRIGHT ! This software was developed within the context of the EUMETSAT Satellite ! Application Facility on Numerical Weather Prediction (NWP SAF), under the ! Cooperation Agreement dated 1 December 2006, between EUMETSAT and the ! Met Office, UK, by one or more partners within the NWP SAF. The partners ! in the NWP SAF are the Met Office, ECMWF, KNMI and MeteoFrance. ! ! Copyright 2010, EUMETSAT, All Rights Reserved. ! ! History: ! Version Date Comment ! ! 1.0 22/07/2010 N.C.Atkinson ! 1.1 21/11/2011 Convert to f90. A. Collard ! ! Code Description: ! FORTRAN 77, following AAPP standards ! ! Declarations: IMPLICIT NONE ! Parameters ! INTEGER(I_KIND), PARAMETER :: nxmax=128 !Max number of spots per scan line ! INTEGER(I_KIND), PARAMETER :: nymax=8192 !Max number of lines. Allows 6hrs of ATMS. INTEGER(I_KIND), PARAMETER :: nxmax=64 !Max number of spots per scan line INTEGER(I_KIND), PARAMETER :: nymax=16384 !Max number of lines. Allows 6hrs of ATMS. REAL(R_KIND) minval, maxval PARAMETER (minval=0.0) !Values less than this are treated as missing PARAMETER (maxval=400.0) !Values greater than this are treated as missing ! Arguments INTEGER(I_KIND), INTENT(IN) :: nx, ny !Size of image REAL(R_KIND), INTENT(INOUT) :: image(nx,ny) !BT or radiance image REAL(R_KIND), INTENT(IN) :: sampling_distx !typically degrees REAL(R_KIND), INTENT(IN) :: sampling_disty !typically degrees REAL(R_KIND), INTENT(IN) :: beamwidth !ditto REAL(R_KIND), INTENT(IN) :: newwidth !ditto REAL(R_KIND), INTENT(IN) :: mtfcutoff !0.0 to 1.0 INTEGER(I_KIND), INTENT(IN) :: nxaverage !Number of samples to average (or zero) INTEGER(I_KIND), INTENT(IN) :: nyaverage !Number of samples to average (or zero) INTEGER(I_KIND), INTENT(IN) :: qc_distx !Number of samples around missing data to set to INTEGER(I_KIND), INTENT(IN) :: qc_disty !Number of samples around missing data to set to INTEGER(I_KIND), INTENT(OUT) :: Error !Error Status ! Local variables INTEGER(I_KIND) :: nxpad, nypad, dx, dy INTEGER(I_KIND) :: i,j,k,ix,iy, ii, jj INTEGER(I_KIND) :: ifirst INTEGER(I_KIND) :: xpow2, ypow2 INTEGER(I_KIND) :: nxav2, nyav2, naverage ! INTEGER(I_KIND) :: deltax INTEGER(I_KIND) :: minii, maxii, minjj, maxjj REAL(R_KIND), ALLOCATABLE :: mtfxin(:),mtfxout(:) REAL(R_KIND), ALLOCATABLE :: mtfyin(:),mtfyout(:) REAL(R_KIND) :: mtfin,mtfout REAL(R_KIND) :: mtfx_constant, mtfy_constant !test REAL(R_KIND), ALLOCATABLE :: mtfpad(:,:) REAL(R_KIND), ALLOCATABLE :: imagepad(:,:) REAL(R_KIND), ALLOCATABLE :: work(:) REAL(R_KIND) :: f,df,factor REAL(R_KIND) :: PI, LN2, LNcsquared LOGICAL :: missing LOGICAL, ALLOCATABLE :: gooddata_map(:,:) ! End of declarations !----------------------------------------- PI = 4.0*atan(1.0) LN2 = LOG(2.0) ! MTF_Constant=-(PI/(2*sampling_dist))**2/LN2 MTFX_Constant=-(PI/(2*sampling_distx))**2/LN2 !test MTFY_Constant=-(PI/(2*sampling_disty))**2/LN2 !test IF (mtfcutoff .GT. 0.0) LNcsquared = LOG(mtfcutoff)**2 nxav2 = nxaverage/2 nyav2 = nyaverage/2 naverage = nxaverage*nyaverage Error = 0 !1) Pad the image up to the nearest power of 2 in each dimension, by reversing !the points near the edge. xpow2 = INT(LOG(nx*1.0)/LN2 + 1.0) ypow2 = INT(LOG(ny*1.0)/LN2 + 1.0) nxpad = 2**xpow2 nypad = 2**ypow2 dx = (nxpad - nx)/2 dy = (nypad - ny)/2 IF (nxpad .GT. nxmax) THEN write(*,*) 'SSMIS_Spatial_Average: nx too large, maximum allowed value is ',nxmax-1 Error = 1 RETURN END IF IF (nypad .GT. nymax) THEN write(*,*) 'SSMIS_Spatial_Average: ny value is ', nypad !test write(*,*) 'SSMIS_Spatial_Average: ny too large, maximum allowed value is ',nymax-1 Error = 1 RETURN END IF ALLOCATE(mtfxin(nxpad),mtfxout(nxpad)) ALLOCATE(mtfyin(nypad),mtfyout(nypad)) ALLOCATE(mtfpad(nxpad,nypad)) ALLOCATE(imagepad(nxpad,nypad)) ALLOCATE(work(nypad)) ALLOCATE(gooddata_map(nxmax,nymax)) !Loop over scan positions DO j=dy+1,dy+ny DO i=dx+1,dx+nx if (image(i-dx,j-dy) > maxval .OR. image(i-dx,j-dy) < minval) & image(i-dx,j-dy) = minval - 1.0_r_kind imagepad(i,j) = image(i-dx,j-dy) !Take a copy of the input data gooddata_map(i,j) = .TRUE. ! Initialised for step 6) !111 format(4i12,2f20.5,l6) ENDDO !Interpolate missing points in the along-track direction ifirst = -1 missing = .false. DO i=dx+1,dx+nx IF (.not.missing) THEN IF (imagepad(i,j) .GE. minval) THEN !imagepad = -1 indicates missing data ifirst = i ELSE missing = .true. ENDIF ELSE IF (imagepad(i,j) .GE. minval) THEN !First good point after missing missing = .false. IF (ifirst .eq. -1) THEN DO k=dx+1,i-1 imagepad(k,j) = imagepad(i,j) !Constant ENDDO ELSE DO k=ifirst+1,i-1 factor = (i-k)*1.0/(i-ifirst) !Interpolate imagepad(k,j) = imagepad(ifirst,j)*factor + & imagepad(i,j)*(1.0-factor) ENDDO ENDIF ENDIF ENDIF ENDDO IF (missing) THEN !Last scan is missing IF (ifirst .GE. 1) then DO k=ifirst+1,dx+nx imagepad(k,j) = imagepad(ifirst,j) !Constant ENDDO ENDIF ENDIF !Continue padding the edges DO i=1,dx imagepad(i,j) = imagepad(dx+dx+2-i,j) ENDDO DO i=nx+dx+1,nxpad imagepad(i,j) = imagepad(nx+dx+nx+dx-i,j) ENDDO ENDDO DO j=1,dy DO i=1,nxpad imagepad(i,j) = imagepad(i,dy+dy+2-j) ENDDO ENDDO DO j=ny+dy+1,nypad DO i=1,nxpad imagepad(i,j) = imagepad(i,ny+dy+ny+dy-j) ENDDO ENDDO !2) Compute the MTF modifications. Assume beams are Gaussian. IF (newwidth .GT. 0) THEN df = 1.0/nxpad DO i=1,nxpad/2+1 f = df*(i-1) !DC to Nyquist ! mtfxin(i) = exp(MTF_Constant*(f*beamwidth)**2) ! mtfxout(i) = exp(MTF_Constant*(f*newwidth)**2) mtfxin(i) = exp(MTFX_Constant*(f*beamwidth)**2) !test mtfxout(i) = exp(MTFX_Constant*(f*newwidth)**2) !test IF (i.GT.1.AND.i.LT.nxpad/2+1) THEN mtfxin(nxpad-i+2) = mtfxin(i) mtfxout(nxpad-i+2) = mtfxout(i) ENDIF ENDDO df = 1.0/nypad DO i=1,nypad/2+1 f = df*(i-1) !DC to Nyquist ! mtfyin(i) = exp(MTF_Constant*(f*beamwidth)**2) ! mtfyout(i) = exp(MTF_Constant*(f*newwidth)**2) mtfyin(i) = exp(MTFY_Constant*(f*beamwidth)**2) !test mtfyout(i) = exp(MTFY_Constant*(f*newwidth)**2) !test IF (i.GT.1.AND.i.LT.nypad/2+1) THEN mtfyin(nypad-i+2) = mtfyin(i) mtfyout(nypad-i+2) = mtfyout(i) ENDIF ENDDO DO i=1,nxpad DO j=1,nypad mtfin = mtfxin(i)*mtfyin(j) mtfout = mtfxout(i)*mtfyout(j) if (mtfcutoff .GT. 0.0) THEN mtfpad(i,j) = (mtfout * & exp(-LN2/LNcsquared*(LOG(mtfout))**2))/mtfin else mtfpad(i,j) = mtfout/mtfin endif ENDDO ENDDO !3) Fourier transform, line by line then column by column. !After each FFT, points 1 to nxpad/2+1 contain the real part of the spectrum, !the rest contain the imaginary part in reverse order. DO j=1,nypad CALL SFFTCF(imagepad(:,j),nxpad,xpow2) ENDDO DO i=1,nxpad DO j=1,nypad work(j) = imagepad(i,j) ENDDO CALL SFFTCF(work,nypad,ypow2) DO j=1,nypad imagepad(i,j) = work(j) ENDDO ENDDO !4) Multiply the spectrum by the MTF factor DO j=1,nypad DO i=1,nxpad imagepad(i,j) = imagepad(i,j)*mtfpad(i,j) ENDDO ENDDO !5) Inverse Fourier transform, column by column then line by line DO i=1,nxpad DO j=1,nypad work(j) = imagepad(i,j) ENDDO CALL SFFTCB(work,nypad,ypow2) DO j=1,nypad imagepad(i,j) = work(j) ENDDO ENDDO DO j=1,nypad CALL SFFTCB(imagepad(:,j),nxpad,xpow2) ENDDO ENDIF !New width is specified !6) Reset missing values in gooddata_map, based on qc_dist and the values ! in the input image array ! Set the ends of the image to missing in the along track direction ! (doing the same across track will remove too much data) ! gooddata_map(:,1:qc_dist)=.FALSE. ! gooddata_map(:,ny-qc_dist+1:ny)=.FALSE. gooddata_map(:,1:qc_disty)=.FALSE. gooddata_map(:,ny-qc_disty+1:ny)=.FALSE. DO j=1,ny DO i=1,nx IF (image(i,j) <= minval) THEN ! minjj=max(j+dy-qc_dist,0) ! maxjj=min(j+dy+qc_dist,nymax) minjj=max(j+dy-qc_disty,0) maxjj=min(j+dy+qc_disty,nymax) DO jj=minjj,maxjj ! deltax=INT(SQRT(REAL(qc_dist**2 - (jj-j-dy)**2 ))) ! minii=max(i+dx-deltax,0) ! maxii=min(i+dx+deltax,nxmax) minii=max(i+dx-qc_distx,0) maxii=min(i+dx+qc_distx,nxmax) DO ii=minii,maxii gooddata_map(ii,jj)=.FALSE. !345 format(10i10) END DO END DO END IF END DO END DO !7) Over-write the input image (points that are not missing) DO j=1,ny DO i=1,nx IF (gooddata_map(i+dx,j+dy)) THEN IF (nxav2 == 0. .AND. nyav2 == 0) THEN image(i,j) = imagepad(i+dx,j+dy) ELSE image(i,j) = 0.0 !Do averaging DO ix = -nxav2,nxav2 DO iy = -nyav2,nyav2 image(i,j) = image(i,j) + imagepad(i+dx+ix,j+dy+iy) ENDDO ENDDO image(i,j) = image(i,j)/naverage ENDIF ELSE image(i,j) = missing_value END IF ENDDO ENDDO !8) Deallocate arrays DEALLOCATE(mtfxin,mtfxout) DEALLOCATE(mtfyin,mtfyout) DEALLOCATE(mtfpad) DEALLOCATE(imagepad) DEALLOCATE(work) DEALLOCATE(gooddata_map) RETURN END SUBROUTINE MODIFY_BEAMWIDTH !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! A real-valued, in place, split-radix FFT program ! Real input and output in data array X ! Length is N = 2 ** M ! Decimation-in-time, cos/sin in second loop ! Output in order: ! [ Re(0), Re(1), ..., Re(N/2), Im(N/2-1), ..., Im(1) ] ! ! This FFT computes ! X(k) = sum_{j=0}^{N-1} x(j)*exp(-2ijk*pi/N) ! ! ! H.V. Sorensen, Rice University, Oct. 1985 ! ! Reference: H.V. Sorensen, D.L. Jones, M.T. Heideman, & C.S. Burrus; ! Real-Valued Fast Fourier Transform Algorithms; IEEE ! Trans. Acoust., Speech, Signal Process.; Vol ASSP-35, ! June 1987, pp. 849-863. ! ! This code was originally named RVFFT. ! ! History: ! 21/11/2011 Converted to something resembling f90. A.Collard ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! SUBROUTINE SFFTCF( X, N, M ) IMPLICIT NONE ! ... Parameters ... REAL(R_KIND), PARAMETER :: SQRT2 = 1.4142135623730950488 REAL(R_KIND), PARAMETER :: TWOPI = 6.2831853071795864769 ! ... Scalar arguments ... INTEGER(I_KIND), INTENT(IN) :: N, M ! ... Array arguments ... REAL(R_KIND), INTENT(INOUT) :: X(N) ! ... Local scalars ... INTEGER(I_KIND) J, I, K, IS, ID, I0, I1, I2, I3, I4, I5, I6, I7, I8 INTEGER(I_KIND) N1, N2, N4, N8 REAL(R_KIND) XT, R1, T1, T2, T3, T4, T5, T6 REAL(R_KIND) A, A3, E, CC1, SS1, CC3, SS3 ! ! ... Exe. statements ... ! IF ( N .EQ. 1 ) RETURN ! J = 1 N1 = N - 1 DO 104 I = 1, N1 IF ( I < J ) THEN XT = X(J) X(J) = X(I) X(I) = XT END IF K = N / 2 DO WHILE ( K < J ) J = J - K K = K / 2 END DO J = J + K 104 CONTINUE ! IS = 1 ID = 4 DO DO 60 I0 = IS, N, ID I1 = I0 + 1 R1 = X(I0) X(I0) = R1 + X(I1) X(I1) = R1 - X(I1) 60 CONTINUE IS = 2 * ID - 1 ID = 4 * ID IF ( IS >= N ) EXIT END DO ! N2 = 2 DO 10 K = 2, M N2 = N2 * 2 N4 = N2 / 4 N8 = N2 / 8 E = TWOPI / N2 IS = 0 ID = N2 * 2 DO DO 38 I = IS, N-1, ID I1 = I + 1 I2 = I1 + N4 I3 = I2 + N4 I4 = I3 + N4 T1 = X(I4) + X(I3) X(I4) = X(I4) - X(I3) X(I3) = X(I1) - T1 X(I1) = X(I1) + T1 IF ( N4 == 1 ) CYCLE I1 = I1 + N8 I2 = I2 + N8 I3 = I3 + N8 I4 = I4 + N8 T1 = ( X(I3) + X(I4) ) / SQRT2 T2 = ( X(I3) - X(I4) ) / SQRT2 X(I4) = X(I2) - T1 X(I3) = - X(I2) - T1 X(I2) = X(I1) - T2 X(I1) = X(I1) + T2 38 CONTINUE IS = 2 * ID - N2 ID = 4 * ID IF ( IS >= N ) EXIT END DO A = E DO 32 J = 2, N8 A3 = 3 * A CC1 = COS(A) SS1 = SIN(A) CC3 = COS(A3) SS3 = SIN(A3) A = J * E IS = 0 ID = 2 * N2 DO DO 30 I = IS, N-1, ID I1 = I + J I2 = I1 + N4 I3 = I2 + N4 I4 = I3 + N4 I5 = I + N4 - J + 2 I6 = I5 + N4 I7 = I6 + N4 I8 = I7 + N4 T1 = X(I3) * CC1 + X(I7) * SS1 T2 = X(I7) * CC1 - X(I3) * SS1 T3 = X(I4) * CC3 + X(I8) * SS3 T4 = X(I8) * CC3 - X(I4) * SS3 T5 = T1 + T3 T6 = T2 + T4 T3 = T1 - T3 T4 = T2 - T4 T2 = X(I6) + T6 X(I3) = T6 - X(I6) X(I8) = T2 T2 = X(I2) - T3 X(I7) = - X(I2) - T3 X(I4) = T2 T1 = X(I1) + T5 X(I6) = X(I1) - T5 X(I1) = T1 T1 = X(I5) + T4 X(I5) = X(I5) - T4 X(I2) = T1 30 CONTINUE IS = 2 * ID - N2 ID = 4 * ID IF ( IS >= N ) EXIT END DO 32 CONTINUE 10 CONTINUE RETURN ! ! ... End of subroutine SFFTCF ... ! END SUBROUTINE SFFTCF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! A real-valued, in place, split-radix IFFT program ! Hermitian symmetric input and real output in array X ! Length is N = 2 ** M ! Decimation-in-frequency, cos/sin in second loop ! Input order: ! [ Re(0), Re(1), ..., Re(N/2), Im(N/2-1), ..., Im(1) ] ! ! This FFT computes ! x(j) = (1/N) * sum_{k=0}^{N-1} X(k)*exp(2ijk*pi/N) ! ! ! H.V. Sorensen, Rice University, Nov. 1985 ! ! Reference: H.V. Sorensen, D.L. Jones, M.T. Heideman, & C.S. Burrus; ! Real-Valued Fast Fourier Transform Algorithms; IEEE ! Trans. Acoust., Speech, Signal Process.; Vol ASSP-35, ! June 1987, pp. 849-863. ! ! This code was originally named IRVFFT. ! ! History: ! 21/11/2011 Converted to something resembling f90. A.Collard ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! SUBROUTINE SFFTCB( X, N, M ) use kinds, only: r_kind,r_double,i_kind IMPLICIT NONE ! ... Parameters ... REAL(R_KIND), PARAMETER :: SQRT2 = 1.4142135623730950488 REAL(R_KIND), PARAMETER :: TWOPI = 6.2831853071795864769 ! ... Scalar arguments ... INTEGER(I_KIND), INTENT(IN) :: N, M ! ... Array arguments ... REAL(R_KIND), INTENT(INOUT) :: X(N) ! ... Local scalars ... INTEGER(I_KIND) J, I, K, IS, ID, I0, I1, I2, I3, I4, I5, I6, I7, I8 INTEGER(I_KIND) N1, N2, N4, N8 REAL(R_KIND) XT, R1, T1, T2, T3, T4, T5 REAL(R_KIND) A, A3, E, CC1, SS1, CC3, SS3 ! ! ... Exe. statements ... ! IF ( N .EQ. 1 ) RETURN ! N2 = 2 * N DO 10 K = 1, M-1 IS = 0 ID = N2 N2 = N2 / 2 N4 = N2 / 4 N8 = N4 / 2 E = TWOPI / N2 DO 17 DO 15 I = IS, N-1, ID I1 = I + 1 I2 = I1 + N4 I3 = I2 + N4 I4 = I3 + N4 T1 = X(I1) - X(I3) X(I1) = X(I1) + X(I3) X(I2) = 2 * X(I2) X(I3) = T1 - 2 * X(I4) X(I4) = T1 + 2 * X(I4) IF ( N4 .EQ. 1 ) CYCLE I1 = I1 + N8 I2 = I2 + N8 I3 = I3 + N8 I4 = I4 + N8 T1 = ( X(I2) - X(I1) ) / SQRT2 T2 = ( X(I4) + X(I3) ) / SQRT2 X(I1) = X(I1) + X(I2) X(I2) = X(I4) - X(I3) X(I3) = 2 * ( - T2 - T1 ) X(I4) = 2 * ( -T2 + T1 ) 15 CONTINUE IS = 2 * ID - N2 ID = 4 * ID IF ( IS >= N-1 ) EXIT END DO A = E DO 20 J = 2, N8 A3 = 3 * A CC1 = COS(A) SS1 = SIN(A) CC3 = COS(A3) SS3 = SIN(A3) A = J * E IS = 0 ID = 2 * N2 DO DO 30 I = IS, N-1, ID I1 = I + J I2 = I1 + N4 I3 = I2 + N4 I4 = I3 + N4 I5 = I + N4 - J + 2 I6 = I5 + N4 I7 = I6 + N4 I8 = I7 + N4 T1 = X(I1) - X(I6) X(I1) = X(I1) + X(I6) T2 = X(I5) - X(I2) X(I5) = X(I2) + X(I5) T3 = X(I8) + X(I3) X(I6) = X(I8) - X(I3) T4 = X(I4) + X(I7) X(I2) = X(I4) - X(I7) T5 = T1 - T4 T1 = T1 + T4 T4 = T2 - T3 T2 = T2 + T3 X(I3) = T5 * CC1 + T4 * SS1 X(I7) = - T4 * CC1 + T5 * SS1 X(I4) = T1 * CC3 - T2 * SS3 X(I8) = T2 * CC3 + T1 * SS3 30 CONTINUE IS = 2 * ID - N2 ID = 4 * ID IF ( IS >= N-1 ) EXIT END DO 20 CONTINUE 10 CONTINUE ! IS = 1 ID = 4 DO DO 60 I0 = IS, N, ID I1 = I0 + 1 R1 = X(I0) X(I0) = R1 + X(I1) X(I1) = R1 - X(I1) 60 CONTINUE IS = 2 * ID - 1 ID = 4 * ID IF ( IS >= N ) EXIT END DO ! J = 1 N1 = N - 1 DO 104 I = 1, N1 IF ( I < J ) THEN XT = X(J) X(J) = X(I) X(I) = XT END IF K = N / 2 DO WHILE ( K < J ) J = J - K K = K / 2 END DO J = J + K 104 CONTINUE XT = 1.0 / FLOAT( N ) DO 99 I = 1, N X(I) = XT * X(I) 99 CONTINUE RETURN ! ! ... End of subroutine SFFTCB ... ! END SUBROUTINE SFFTCB END MODULE SSMIS_Spatial_Average_Mod