program create_daily_global_csrbs c Code to read HIRS orbital CSRB files and create global, 1.0-degree gridded c statistics. c------------------------------------------------------------------------------- implicit none c------------------------------------------------------------------------------- integer nlon_bins parameter (nlon_bins = 360) integer nlat_bins parameter (nlat_bins = 180) integer nstats parameter (nstats = 4) integer nbands parameter (nbands = 8) integer nscns parameter (nscns = 3) integer reclen parameter (reclen = nlon_bins * nlat_bins * nbands * nscns * nstats * 8) integer nb, nt parameter (nb=20,nt=2) real cwn(nb), fk1(nb), fk2(nb), tc(nt,nb) real cwn_shifted(nb), fk1_shifted(nb), fk2_shifted(nb), tc_shifted(nt,nb) common/hirpfc/cwn, fk1, fk2, tc common/hirpfc_shifted/cwn_shifted, fk1_shifted, fk2_shifted, tc_shifted c------------------------------------------------------------------------------- character*120 Input_File, Output_File, parm character*6 sat integer i, j, k, m, ierr, mbnds(nbands), shifted_FM_opt real*8 n(nlon_bins, nlat_bins, nbands, nscns) real*8 sum_obsclr(nlon_bins, nlat_bins, nbands, nscns) real*8 sum_calclr(nlon_bins, nlat_bins, nbands, nscns) real*8 ssq_raddiff(nlon_bins, nlat_bins, nbands, nscns) real*8 mean_obsclr real*8 mean_obsclr_bt real*8 mean_calclr real*8 mean_calclr_bt real*8 sum_diff, sqsum_diff, num, den real*4 mean_diff(nlon_bins, nlat_bins, nbands, nscns) real*4 mean_diff_bt(nlon_bins, nlat_bins, nbands, nscns) real*4 sdev(nlon_bins, nlat_bins, nbands, nscns) real*4 nout(nlon_bins, nlat_bins, nbands, nscns) real*4 radobs, radcal c------------------------------------------------------------------------------- c External functions. real rf_hirsbright external rf_hirsbright, rf_hirspfco_101 c------------------------------------------------------------------------------- c Data statements. data mbnds /4,5,6,7,8,10,11,12/ c------------------------------------------------------------------------------- c Get input file name. call getarg(1,parm) read(parm,'(a120)') Input_File write(*,'(2x,'' Input File is: '',A100)') Input_File c------------------------------------------------------------------------------- c Get output file name. call getarg(2,parm) read(parm,'(a120)') Output_File write(*,'(2x,'' Output File is: '',A100)') Output_File c------------------------------------------------------------------------------- c Get input satellite number. call getarg(3,parm) read(parm,'(a6)') sat write(*,'(2x,'' Satellite #: '',a6)') sat c------------------------------------------------------------------------------- c Get input shifted forward model option. call getarg(4,parm) read(parm,'(i4)') shifted_FM_opt write(*,'(2x,'' Shifted FM option: '', i5)') shifted_FM_opt c------------------------------------------------------------------------------- c Open and read input file. open(30,file=Input_File,form='unformatted',access='direct',iostat=ierr, + recl=reclen) if(ierr .ne. 0) then write(*,'(2x,''Could not open input data file '', + A70,i10)') Input_File, ierr go to 100 end if read(30, rec=1, iostat=ierr) n, sum_obsclr, sum_calclr, ssq_raddiff if(ierr .ne. 0) then write(*,'(2x,''Could not read input data file '', + A70,i10)') Input_File, ierr close(30) go to 100 end if close(30) c------------------------------------------------------------------------------- c Get HIRS instrument data for current satellite. call rf_hirspfco_101(sat, shifted_FM_opt) c------------------------------------------------------------------------------- c Loop through the data. do m = 1, nscns do k = 1, nbands do j = 1, nlat_bins do i = 1, nlon_bins c------------------------------------------------------------------------------- c Calculate daily means and differences from input sums. nout(i,j,k,m) = real(n(i,j,k,m)) if(n(i,j,k,m) .gt. 0.d0) then mean_obsclr = sum_obsclr(i,j,k,m) / n(i,j,k,m) mean_calclr = sum_calclr(i,j,k,m) / n(i,j,k,m) mean_diff(i,j,k,m) = real(mean_obsclr - mean_calclr) radobs = real(mean_obsclr) mean_obsclr_bt = rf_hirsbright(radobs, mbnds(k), 1) radcal = real(mean_calclr) mean_calclr_bt = rf_hirsbright(radcal, mbnds(k), shifted_FM_opt) mean_diff_bt(i,j,k,m) = real(mean_obsclr_bt - mean_calclr_bt) if(mean_diff(i,j,k,m) .gt. 20.0) then write(*,'(4i4,10f10.3)') i,j,k,mbnds(k),n(i,j,k,m), * sum_obsclr(i,j,k,m),mean_obsclr,mean_calclr,mean_diff(i,j,k,m), * radobs,mean_obsclr_bt,radcal,mean_calclr_bt, * mean_diff_bt(i,j,k,m) end if else mean_diff(i,j,k,m) = -999.0 mean_diff_bt(i,j,k,m) = -999.0 end if c------------------------------------------------------------------------------- c Calculate daily standard deviation of differences. if(n(i,j,k,m) .gt. 1.d0) then sum_diff = n(i,j,k,m) * mean_diff(i,j,k,m) sqsum_diff = sum_diff * sum_diff num = (n(i,j,k,m)*ssq_raddiff(i,j,k,m)) - sqsum_diff den = n(i,j,k,m) * (n(i,j,k,m) - 1) if(num .gt. 0.d0) then sdev(i,j,k,m) = real(dsqrt(num / den)) else sdev(i,j,k,m) = 0.0 end if else if(n(i,j,k,m) .eq. 1.d0) then sdev(i,j,k,m) = 0.0 else sdev(i,j,k,m) = -999.0 end if enddo enddo enddo enddo c------------------------------------------------------------------------------- c Output stats to binary file. open(40, file=Output_File, form='unformatted', * access='sequential') write(40) nout, mean_diff, mean_diff_bt, sdev close(40) c------------------------------------------------------------------------------- 100 continue end