program ukmhires !----------------------------------------------------------------- ! Merges ukm hires data in 8 regional patches into one global array ! Data Source: /dcom/us007003/yyyymmdd/wgrbbul/ukmet_hires ! Use bi-linear interpolation to create a 1-deg global dataset. ! Fanglin Yang, June 2014 !----------------------------------------------------------------- !-------------------------------------- !---input data info !-------------------------------------- integer, parameter :: np=8 !number of patches integer, parameter :: nx=108 !number of longtudinal points for each patch integer, parameter :: ny=162 !number of latitudinal points for each patch integer, parameter :: nxny=nx*ny !total points for each patch integer, parameter :: mmax=200000 !maximum number of points to unpack real, parameter :: resy=556.0, resx=833.0 !data resolution (x1000) logical*1 :: lb(mmax) character*200 :: input,output integer :: iargc external :: iargc integer :: nargs ! number of command-line arguments character*200 :: argument ! space for command-line argument integer :: jpds(200),jgds(200), fhr integer :: kpds(200),kgds(200) real :: varp(mmax,np) !data from each patch real :: var(nx*4+1,ny*2) !a global high-res array real :: xs(np), xe(np), ys(np), ye(np) !starting and ending lat-lon of patches real :: lon(nx*4+1),lat(ny*2) !lat-lon positions !-------------------------------------- !---output data, 1x1 deg global !-------------------------------------- real, parameter :: resg=1000.0 integer, parameter :: mx=360000/resg,my=180000/resg+1 integer, parameter :: mxmy=mx*my real :: varg(mx,my),varg1(mx*my) real :: long(mx),latg(my) real :: dy1(my),dy2(my),dx1(mx),dx2(mx) integer :: ix(mx),iy(my) real :: tmp1,tmp2,sums,sumn !------------------------------------------------------------------------------------- !--starting and ending lat and lon of each patch (x1000) data xs/341250, 71250, 161250, 251250, 341250, 71250, 161250, 251250/ data xe/ 70416, 160416, 250416, 340416, 70416, 160416, 250416, 340416/ data ys/ 279, 279, 279, 279, -89721, -89721, -89721, -89721/ data ye/ 89722, 89722, 89722, 89722, -278, -278, -278, -278/ !---source data position do n=1,4 i0=(n-1)*nx do i=i0+1,i0+nx lon(i)=xs(n)+resx*(i-i0-1) if(lon(i).gt.360000) lon(i)=lon(i)-360000 enddo !!print*,n,xe(n),lon(i0+nx) enddo lon(4*nx+1)=lon(1) do j=1,ny lat(j)=ys(5)+resy*(j-1) !southern hemisphere enddo do j=ny+1,ny*2 lat(j)=ys(1)+resy*(j-ny-1) !northern hemisphere enddo !---output data position do i=1,mx long(i)=(i-1)*resg enddo do j=1,my latg(j)=-90000.0+(j-1)*resg enddo !---temporary arrays used for linear interpolation do j=1,my do jj=1,2*ny-1 if(latg(j).ge.lat(jj).and.latg(j).lt.lat(jj+1)) then dy1(j)=latg(j)-lat(jj) dy2(j)=lat(jj+1)-latg(j) iy(j)=jj endif enddo enddo !special case at 0E i=1; ii=23; ix(i)=23 dx1(i)=360000.0-lon(ii) dx2(i)=lon(ii+1) do i=2,mx do ii=1,4*nx if(long(i).ge.lon(ii).and.long(i).lt.lon(ii+1)) then dx1(i)=long(i)-lon(ii) dx2(i)=lon(ii+1)-long(i) ix(i)=ii endif enddo enddo !----------------------------------------------------- nargs = iargc() ! iargc() - number of arguments if (nargs.lt.3) then write(*,*)'usage : ukm_hires_merge.x input output fcst_hour' stop endif call getarg(1,argument) read(argument,*) input call getarg(2,argument) read(argument,*) output call getarg(3,argument) read(argument,*) fhr print*, trim(input)," ",trim(output), " fcst hour=",fhr fhrm6=fhr-6 fhrm12=fhr-12 call baopenr(10,trim(input),iret) if (iret .ne. 0) write(6,*)" failed to open ", input call baopen(20,trim(output),iret) if (iret .ne. 0) write(6,*)" failed to open ", output if (iret .ne. 0) goto 200 nrec=-1 10 continue jpds=-1; jgds=-1 do n=1,np jgds(2)=nx; jgds(3)=ny; jgds(9)=resx; jgds(10)=resy jgds(4)=ys(n); jgds(5)=xs(n); jgds(7)=ye(n); jgds(8)=xe(n) call getgb(10,0,mmax,nrec,jpds,jgds,kf,k,kpds,kgds,lb,varp(1,n),iret) if(iret.ne.0) goto 100 !reached end of record or incorrect file ! write(1,*) ! write(1,'("patch=",i3," iret=",i2," irec=",i10," nxny=",i10)') n,iret,k,kf ! write(1,'("jpds=",16i8)') (jpds(k),k=1,16) ! write(1,'("kpds=",16i8)') (kpds(k),k=1,16) ! write(1,'("jgds=",16i8)') (jgds(k),k=1,16) ! write(1,'("kgds=",16i8)') (kgds(k),k=1,16) !--look for the same variable in different patches and for the same forecast time do k=1,16 jpds(k)=kpds(k) enddo nrec=nrec-1 enddo !---assemble a global arrary for right forecast hour !-- for precipitation APCP, jpds(15)=fhr if((jpds(14).eq.fhr) .or. (jpds(5).eq.61 .and. jpds(15).eq.fhr) ) then do n=1,4 !patches in northern hemisphere do j=1,ny do i=1,nx jj=ny+j ii=(n-1)*nx+i var(ii,jj)=varp((j-1)*nx+i,n) enddo enddo enddo do n=5,8 !patches in southern hemisphere do j=1,ny do i=1,nx jj=j ii=(n-5)*nx+i var(ii,jj)=varp((j-1)*nx+i,n) enddo enddo enddo do jj=1,2*ny var(4*nx+1,jj)=var(1,jj) enddo !--interpolate to resg-deg global array !--for north pole and south pole sums=0.0; sumn=0.0 do i=1,4*nx sums=sums+var(i,1) sumn=sumn+var(i,2*ny) enddo do i=1,mx varg(i,1)=sums/(4*nx) varg(i,my)=sumn/(4*nx) enddo !--for other points do j=2,my-1 do i=1,mx tmp1=(var(ix(i),iy(j))*dy2(j)+var(ix(i),iy(j)+1)*dy1(j))/(dy1(j)+dy2(j)) tmp2=(var(ix(i)+1,iy(j))*dy2(j)+var(ix(i)+1,iy(j)+1)*dy1(j))/(dy1(j)+dy2(j)) varg(i,j)=(tmp1*dx2(i)+tmp2*dx1(i))/(dx1(i)+dx2(i)) enddo enddo do j=1,my do i=1,mx varg1((j-1)*mx+i)=varg(i,j) enddo enddo !--write out grib1 global array kgds(2)=mx; kgds(3)=my; kgds(9)=resg; kgds(10)=resg kgds(4)=-90000; kgds(5)=0; kgds(7)=90000; kgds(8)=(mx-1)*resg call putgb(20,mxmy,kpds,kgds,lb,varg1,iret) endif goto 10 100 call baclose(10,iret) call baclose(20,iret) 200 end