program cal_swath_binary !!! read hourly or 0.5hrly output of precip and wind10m from NMMB implicit none integer :: i,j,l, k character(len=256) :: infile,outfile integer :: im2,jm2,lm2, km2 real, allocatable, dimension(:,:) :: totalrain2 real, allocatable, dimension(:,:) :: hlat2,hlon2 real, allocatable, dimension(:,:) :: u2,v2,w2 ! 10-m wind speed integer :: im3,jm3,lm3, km3 real, allocatable, dimension(:,:) :: totalrain3 real, allocatable, dimension(:,:) :: hlat3,hlon3 real, allocatable, dimension(:,:) :: u3,v3,w3 ! 10-m wind speed integer :: is,js real, allocatable, dimension(:,:) :: rain,wind10 ! rain and wind10 , swath real, allocatable, dimension(:) :: xx,yy ! lat, lon for swath domain integer :: itime, itmax_track, iloop, kmax2,kmax3 real :: dt,tmax, tt, tmpdis, tmns real fhr(50),rlat(50),rlon(50) ! track information real :: tmpwind,tmppress, x0,x1,y0,y1,dxy,stormsize real :: dlat,dlon parameter(dxy=0.03,stormsize=6.0 ) !parameter(dxy=0.2,stormsize=6.0) integer :: cross180, i8,j8 ! ! read short track file, fort.21 ! do i=1,50 read(21,*,end=199)fhr(i),rlat(i),rlon(i),tmpwind,tmppress print*,'hr,lat,lon=',fhr(i),rlat(i),rlon(i) enddo 199 itmax_track=i-1 ! determin swath domain ! dxy=0.05 ! resolution ! corners y0=minval(rlat(1:itmax_track))-stormsize y1=maxval(rlat(1:itmax_track))+stormsize x0=minval(rlon(1:itmax_track))-stormsize x1=maxval(rlon(1:itmax_track))+stormsize print*,'min,max LON=',x0,x1 print*,'min,max LAT=',y0,y1 !! Weiguo Wang , fix for data cross dateline cross180=0 if ( abs(x1-x0) > 180.0 .and. x0 < 0 ) then cross180=1 do i=1,itmax_track if(rlon(i) > 0)rlon(i)=rlon(i)-360.0 enddo x0=minval(rlon(1:itmax_track))-stormsize x1=maxval(rlon(1:itmax_track))+stormsize print*,'found track cross dateline' print*,'min,max LON=',x0,x1 endif is=nint( (x1-x0)/dxy ) js=nint( (y1-y0)/dxy ) print*,'is,js=',is,js ! swath array allocate (wind10(is,js),rain(is,js),yy(js),xx(is) ) do i=1,is xx(i)=(i-1)*dxy + x0 enddo do j=1,js yy(j)=(j-1)*dxy + y0 enddo ! initialize array wind10=-999.0 rain=-999.0 ! distance do k=1,itmax_track do i=1,is do j=1,js tmpdis=( xx(i)-rlon(k) )**2 + (yy(j)-rlat(k))**2 if( tmpdis .le. stormsize*stormsize ) then wind10(i,j)=0.0 rain(i,j)=0.0 endif enddo enddo enddo ! !!!!!!!!!!!!!!! get dim of saved binary data !!!!!!!!!!!!! call rewrite('swath_d02','swath_d02r',im2,jm2,kmax2) call rewrite('swath_d03','swath_d03r',im3,jm3,kmax3) !! estimate save interval,dt, from kmax, total integration is 126 hr. dt = 126.0 / (kmax2-1) open(22,file='swath_d02r',form='unformatted',CONVERT="BIG_ENDIAN") open(23,file='swath_d03r',form='unformatted',CONVERT="BIG_ENDIAN") !! allocate variables allocate (hlat2(im2,jm2),hlon2(im2,jm2)) allocate (totalrain2(im2,jm2)) allocate (w2(im2,jm2)) allocate (hlat3(im3,jm3),hlon3(im3,jm3)) allocate (totalrain3(im3,jm3)) allocate (w3(im3,jm3)) !do i=1,kmax2 ! kmax2=kmax3 do i = 1, min(kmax2, ifix( fhr(itmax_track)/dt+1 ) ) ! kmax2=kmax3 !D2 print*,'i=,time=',i,dt*(i-1),fhr(itmax_track) print*,'reading...swath_d02r' read(22)hlat2 read(22)hlon2 read(22)totalrain2 read(22)w2 hlat2 = hlat2*180.0/3.14159265359 hlon2 = hlon2*180.0/3.14159265359 !wang if (cross180 > 0 ) then do i8=1,im2 do j8=1,jm2 if(hlon2(i8,j8) > 0) hlon2(i8,j8) = hlon2(i8,j8) -360.0 enddo enddo endif ! !D3 print*,'reading...swath_d03r' read(23)hlat3 read(23)hlon3 read(23)totalrain3 read(23)w3 hlat3 = hlat3*180.0/3.14159265359 hlon3 = hlon3*180.0/3.14159265359 !wang if (cross180 > 0 ) then do i8=1,im3 do j8=1,jm3 if(hlon3(i8,j8) > 0) hlon3(i8,j8) = hlon3(i8,j8) -360.0 enddo enddo endif ! ! print*,'min,max hlon1',minval(hlon1),maxval(hlon1) ! print*,'min,max hlat1',minval(hlat1),maxval(hlat1) ! call interp_homo(hlon1,hlat1,im1,jm1,w1 & ! ,xx,yy,is,js,wind10 ) call interp_homo(hlon2,hlat2,im2,jm2,w2 & ,xx,yy,is,js,wind10 ) call interp_homo(hlon2,hlat2,im2,jm2,totalrain2 & ,xx,yy,is,js,rain ) if (1 == 2 .and. i > 1 .and. i < itmax_track) then ! interpolate hourly tt=0.5 ! hourly dlat=rlat(i)-rlat(i-1) dlon=rlon(i)-rlon(i-1) do iloop=1,dt/tt -1 ! move dom (dt/tt-1) times hlon2=hlon2+dlon/dt * tt hlat2=hlat2+dlat/dt * tt call interp_homo(hlon2,hlat2,im2,jm2,w2 & ,xx,yy,is,js,wind10 ) call interp_homo(hlon2,hlat2,im2,jm2,totalrain2 & ,xx,yy,is,js,rain ) enddo endif call interp_homo(hlon3,hlat3,im3,jm3,w3 & ,xx,yy,is,js,wind10 ) call interp_homo(hlon3,hlat3,im3,jm3,totalrain3 & ,xx,yy,is,js,rain ) enddo ! time close(22) close(23) ! write(*,*)'im,jm=',im1,jm1 ! write(*,*)'xx=',xx(1:10) ! write(*,*)'yyy=',yy(1:10) ! write(1,*)totalrain1 ! write(3,*)w1 ! write(4,*)wind10 ! write(*,*)'nint(3.6),nint(3.4),nint(-3.4)nint(-3.6)=',nint(3.6),nint(3.4),nint(-3.4),nint(-3.6) ! write(*,*)'ifix(3.6),ifix(3.4),ifix(-3.4),ifix(-3.6)=',ifix(3.6),ifix(3.4),ifix(-3.4),ifix(-3.6) !!!! OUTPUT in GRaDs format open(80,file='swath.dat', & form="unformatted") !,access='direct',recl=nx*ny) ! write(80)crain_obs,crain_mod,width_obs,width_mod,cx,cy write(80)wind10, rain close(80) !! write(82,*)wind10, rain open(81,file='swath.ctl') write(81,*)"dset ^swath.dat" write(81,*)"undef ",-999.0 write(81,*)"xdef ",is," linear ",xx(1),dxy write(81,*)"ydef ",js," linear ",yy(1),dxy write(81,*)"zdef 1 LEVELS 0" write(81,*)"tdef 1 linear 00Z01JAN2011 1mon" write(81,*)"vars 2" write(81,*)"wind10 0 99 10-m wind m/s" write(81,*)"rain 0 99 rain m" write(81,*)"endvars" close(81) ! now write data out to text file for NHC ! make the ascii files for Michelle's color plots write(82,3332)xx(1),xx(is),yy(1),yy(js),dxy,is,js write(83,3332)xx(1),xx(is),yy(1),yy(js),dxy,is,js 3332 format(5f9.2,2i5) do j = 1 , js do i = 1 , is write(82,3333) yy(j),xx(i),wind10(i,j) write(83,3333) yy(j),xx(i),rain(i,j) 3333 format(3f9.2) enddo enddo tmns = -99.0 write(82,3333)tmns,tmns,tmns write(83,3333)tmns,tmns,tmns ! include duration & track of storm write(82,3335)ifix( fhr(itmax_track) ) write(83,3335)ifix( fhr(itmax_track) ) 3335 format(i5) do i = 1 , itmax_track write(82,3334) rlat(i),rlon(i) write(83,3334) rlat(i),rlon(i) 3334 format(2f8.2) enddo end program cal_swath_binary subroutine interp_homo(xold,yold,nxold,nyold,indata & ,xnew,ynew,nxnew,nynew,outdata ) ! interpolate data from ununiform grid to uniform grid ! implicit none integer :: nxold,nyold,nxnew,nynew real,dimension(nxold,nyold) :: xold,yold,indata real,dimension(nxnew) :: xnew real,dimension(nynew) :: ynew real,dimension(nxnew,nynew) :: outdata !local integer i,j,k ,ii,jj ,ki,kj real,dimension(nxnew,nynew) :: tmp1,tmp2 real dx, dy, sum1,sum2, dis2, dis2old, x1,x2, y1,y2 dx=xnew(2)-xnew(1) dy=ynew(2)-ynew(1) x1=minval(xold) x2=maxval(xold) y1=minval(yold) y2=maxval(yold) !!! outdata(:,:)=-999.0 tmp1(:,:)=0.0 tmp2(:,:)=0.0 do i=1,nxold do j=1,nyold !position of old x/y in new system ii=nint( (xold(i,j) - xnew(1))/dx ) + 1 jj=nint( (yold(i,j) - ynew(1))/dy ) + 1 if( ii .ge. 1 .and. ii .le. nxnew .and. & jj .ge. 1 .and. jj .le. nynew ) then dis2= ( xold(i,j)-xnew(ii) )**2 + ( yold(i,j) - ynew(jj) )**2 if (dis2 < 1e-5) dis2=1.e-5 ! dis2=1.0 tmp1(ii,jj)=tmp1(ii,jj)+1.0/dis2 tmp2(ii,jj)=tmp2(ii,jj)+indata(i,j)/dis2 endif enddo enddo do i=1,nxnew do j=1,nynew if(tmp1(i,j) > 0.0) then tmp2(i,j)=tmp2(i,j)/tmp1(i,j) else tmp2(i,j)=-999.0 endif enddo enddo ! fill -999, when old grid is coarser than new, some points may not have values ! tmp2 save interpolated values at this call tmp1=tmp2 ! save out data do i=1,nxnew do j=1,nynew if (tmp1(i,j) == -999. .and. & ! make sure interpolated points are within old domain. No extroplation xnew(i) .ge. x1 .and. xnew(i) .le. x2 .and. & ynew(j) .ge. y1 .and. ynew(j) .le. y2 ) then sum1=0.0 sum2=0.0 dis2old=9.0e20 do ki=-4,4 do kj=-4,4 ii=i+ki jj=j+kj if (ii .ge.1 .and. ii .le. nxnew .and. & jj .ge.1 .and. jj .le. nynew .and. & tmp1(ii,jj) /= -999.0 ) then dis2=(xnew(i)-xnew(ii))**2+(ynew(j)-ynew(jj))**2 if (dis2 < 1.e-5) dis2=1.0e-5 sum1=sum1+ 1.0/dis2 sum2=sum2+ tmp1(ii,jj)/dis2 ! find nearest point ! if(dis2 < dis2old)then ! dis2old=dis2 ! sum1=1.0 ! sum2=tmp1(ii,jj) ! endif ! find nearst point endif enddo !ki enddo ! kj if(sum1 /=0.0) tmp2(i,j)=sum2/sum1 endif enddo !j enddo !i ! to build swath plot, compare previous values (outdata) with tmp2, then output do i=1,nxnew do j=1,nynew if(tmp2(i,j) /= -999.0 & .and. tmp2(i,j) >= outdata(i,j) .and. outdata(i,j) /= -999.0) then outdata(i,j)=tmp2(i,j) endif enddo enddo end subroutine interp_homo subroutine rewrite(infile,outfile,nx,ny,kmax) ! rewrite data so that they can be read hour by hour character(len=*)infile, outfile integer nx,ny,kmax, k real*4, dimension(:,:,:), allocatable::glat,glon,prec,w10m open(10,file=trim(infile),form='unformatted',convert='big_endian',status='old') read(10)nx,ny,kmax print*,'nx,ny,kmax=',nx,ny,kmax allocate(glat(nx,ny,kmax)) allocate(glon(nx,ny,kmax)) allocate(prec(nx,ny,kmax)) allocate(w10m(nx,ny,kmax)) do k=1,kmax read(10)glat(:,:,k) enddo do k=1,kmax read(10)glon(:,:,k) enddo do k=1,kmax read(10)prec(:,:,k) print*,'min,max prec=',minval(prec(:,:,k)),maxval(prec(:,:,k)) enddo do k=1,kmax read(10)w10m(:,:,k) print*,'min,max wind10m=',minval(w10m(:,:,k)),maxval(w10m(:,:,k)) enddo close(10) open(20,file=trim(outfile),form='unformatted',convert='big_endian') do k=1,kmax print*,'writing k=',k write(20)glat(:,:,k) write(20)glon(:,:,k) write(20)prec(:,:,k) write(20)w10m(:,:,k) enddo close(20) deallocate(glat) deallocate(glon) deallocate(prec) deallocate(w10m) end subroutine rewrite