PROGRAM RMORPH c c original code at: /export-3/cpclnxsvr1/cpcsat/rmorph/version1 c c c c Adjusts CMORPH half-hourly, 8km precip. by the ratio of daily US c daily gauge precip & daily CMORPH (12Z-12Z) which is c functionally the same as disaggregating, i.e.: c c CMORPH half-hourly c ------------------ X Gauge is equivalent to: c CMORPH daily c c c Gauge c . CMORPH half-hourly X ----- c CMORPH daily c c c c c It's important that the c CMORPH data match the daily accumulation period. For c example, the US & S. American gauge data are 12Z-12Z accumulations. c c This version uses a new adjustment scheme over South America. c Instead of disaggregating the daily gauge data, this new c scheme multiplies each cmorph value by the ratio of the daily 1x1 S. American c gauge value and the 24-h cmorph accumulation over that same 1x1 grid c domain. NOTE: this scheme uses 9-pt smoothed values of the c ratios c c By doing it this way, the adjusted cmorph values ("rmorph") c retain the spatial variability within the 1x1 deg grid, but when c the individual values within the 1x1 deg. gridbox are summed over the c daily period, the sum matches the gauge analysis value. View c test results in "test/test_method.f" c c NOTE: Program (particularly subroutine wtrmorph) is set up for 12z-12z c periods, so modifications will be necessary for daily rain gauge data sets that c use different 24-h periods for reporting daily rainfall totals. c c-------------------------------------------------------------------------------- c c PARAMETER (ius=321, jus=201, isa=61,jsa=76, icx=4948, jcy=1649) real*4 gag(ius,jus), cm(icx,jcy,48), xptsa(isa,jsa), x usmask(ius,jus), tempus(ius,jus), gsa(isa,jsa), x samask(isa,jsa), tempsa(isa,jsa), cmsa(isa,jsa), x cmus(ius,jus), ratio_sa(isa,jsa), xptsus(ius,jus), x ratio_smth_sa(isa,jsa), gnum(isa,jsa), x ratio_us(ius,jus), ratio_smth_us(ius,jus) character*1 buf(icx,jcy) character*120 f_gauge data io_gauge /1/, io_mask /2/, iosag /50/ open (77,file='diag.us.mask',recl=ius*jus*4,access='direct') open (78,file='diag.us.data',recl=ius*jus*4,access='direct') open (97,file='diag.sa.mask',recl=isa*jsa*4,access='direct') open (98,file='diag.sa.data',recl=isa*jsa*4,access='direct') open (99,file='diag.print') c c................................................................ c get land/sea masks for SA & US call masksa(io_mask,isa,jsa,tempsa,samask) write(97,rec=1) ((samask(i,j),i=1,isa),j=1,jsa) call maskus(io_mask,ius,jus,usmask) write(77,rec=1) ((usmask(i,j),i=1,ius),j=1,jus) c................................................................ c read(*,*) iyrlast,imolast,idalast,iyr,imo,nday1,nday2 idatlast=iyrlast*10000+imolast*100+idalast nday=0 do 500 ida=nday1,nday2 iyrmndy=iyr*10000+imo*100+ida nday=nday+1 c c...... Read in daily US gauge analysis data call gauge_us(io_gauge,iyrmndy,ius,jus,gag, x usmask,tempus,iretus) write(78,rec=(nday-1)*4+1) gag c c c c...... Read in daily S. America gauge analysis data; also, read c raingauge station file & create array ("gnum") that contians the # of c rain gauge reports for each gridbox in the S. America daily rain gauge c analysis c grdint=1.0 ntot=0 do j=1,ny do i=1,nx gnum(i,j)=0. enddo enddo call gauge_sa(io_gauge,iyrmndy,isa,jsa,gsa, x samask,tempsa,iretsa) call gtname(iyr,imo,ida,f_gauge,iret) c if(iret.eq.1) then call grid(iosag,f_gauge,grdint,isa,jsa,ntot,gnum) else write(99,*) 'No SA gauge file for: ',iyr,imo,ida endif c write(98,rec=(ida-1)*6+1) gsa write(98,rec=(ida-1)*6+2) gnum c cNOTE ....... N & S America gauge analyses are 12Z-12Z c jhr=0 call gtcmorph(iyrlast,imolast,idalast,12,23,icx,jcy,buf,jhr,cm) jhr=12 call gtcmorph(iyr,imo,ida,00,11,icx,jcy,buf,jhr,cm) c c...... Work on gauge adjustment over US c if(iretus.eq.0) then xlonw=219.9070 c cc ylatn=59.963614 c....... The following value is 59.963614 + 0.125 ... 0.125 is added to center c the 1/4 deg data w.r.t. 8km CMORPH data so that the gauge c adjustment is applied appropriately ylatn=60.451114 c call mk_ratio_us (ylatn,xlonw,cm,icx,jcy,usmask,gag,cmus, x xptsus,ratio_us,ius,jus,nday) c thresh=7.0 call adjust_us (cm,icx,jcy,usmask,xlonw,ylatn,gag,cmus, x ratio_us,ius,jus,thresh,ida) write(99,*) 'Thru with US gauge adjustment' else write(99,*)'No US gauge data for this date' endif c c...... Work on gauge adjustment over South America c if(iretsa.eq.0) then xlonw=270.125 ylatn=14.875 call mk_ratio_sa (ylatn,xlonw,cm,icx,jcy,samask,gsa,cmsa, x xptsa,ratio_sa,isa,jsa,nday) call smth9(ratio_sa,gnum,isa,jsa,ratio_smth_sa) write(98,rec=(ida-1)*6+6) ratio_smth_sa thresh=7.0 call adjust_sa (cm,icx,jcy,samask,xlonw,ylatn,gsa,cmsa, x ratio_smth_sa,gnum,isa,jsa,thresh,ida) write(99,*) 'Thru with SA gauge adjustment' else write(99,*)'No SA gauge data for this date' endif call wtrmorph(iyrlast,imolast,idalast,iyr,imo,ida,icx,jcy,buf,cm) iyrlast=iyr imolast=imo idalast=ida write(99,*) '------------- done with date:',iyrmndy write(*,*) '------------- done with date:',iyrmndy 500 continue stop end c c---------------------------------------------------------------- c subroutine adjust_sa (cm,icx,jcy,samask,xlonw,ylatn,gsa,cmsa, x ratio,gnum,isa,jsa,thresh,ida) real*4 cm(icx,jcy,48), gsa(isa,jsa), cmsa(isa,jsa), x samask(isa,jsa), ratio(isa,jsa), gnum(isa,jsa), h(48) c c c "ii" & "jj" are indices that map the Cmorph data to the c appropriate S. American analysis gridbox. c i1=3713 ! 270.125 i2=4538 ! 330.125 c j1= 621 ! 14.875N Note: CMORPH oriented N-S so enter "-14.875" in Grads to get this value j2=1649 ! 59.875S c c c..... Adjust the half-hourly, 8km CMORPH values by the ratio between c the 1x1 deg gauge & 1x1 deg CMORPH c xbeg=0.036378335 xint=0.072756669 c ybeg=59.963614 yint=0.072771377 do 300 j=j1,j2 ylat=ybeg - (j-1)*yint jj=ylatn-ylat+1. if(jj.lt.1) jj=1 if(jj.gt.jsa) jj=jsa do 200 i=i1,i2 xlon=(i-1)*xint + xbeg ii=xlon-xlonw + 1 if(ii.lt.1) ii=1 if(ii.gt.isa) ii=isa c if(xlon.ge.310.7 .and. xlon.le.310.8 .and. ylat.ge.11.7.and. c x ylat.le.11.8) then if(samask(ii,jj).lt.1.) then ccccc write(99,*) 'SA mask < 1 at: ',xlon,ylat,ii,jj go to 200 endif if(gsa(ii,jj).lt.0.) go to 200 do 100 l=1,48 h(l)=cm(i,j,l) c c.............. If daily CMORPH total is nonzero, then, if the hourly c CMORPH is > zero, multiply each hourly CMORPH value by the ratio c between the daily gauge and daily CMORPH totals c if(cmsa(ii,jj).gt.0.) then if(cm(i,j,l).gt.0.) then cm(i,j,l)=cm(i,j,l) * ratio(ii,jj) endif go to 90 endif c c.............. If daily CMORPH total is < or = 0 (i.e. all hourly values are c missing or zero, and the daily gauge total > "thresh", set all c hourly CMORPH values to 1/24th of the daily gauge amt.; otherwise, set c all CMORPH hourly amts. to zero -- this is done because the rain gauge c amounts are small and may be spurious due to the rain c gauge objective analysis technique, i.e. spreading out of low precip c amounts due to smoothing. c if(cmsa(ii,jj).le.0.) then if(gnum(ii,jj).eq.0.) go to 200 if(gsa(ii,jj).gt.thresh) then c c................... Because the final RMORPH values are scaled (by "5") c and stored as a scaled integer in a byte, values of c hourly precip that are < 1 will be evaluated as zero. Therefore, if the c daily rain gauge value is < 4.8 (i.e. 24/5), then instead of setting all 24 c hourly values to 1/24th of the daily gauge value for this situation, set c the values of the 4 synoptic times to 1/4th of the daily gauge value. c if(gsa(ii,jj) .lt. 4.8) then cm(i,j,1)= 6.*gsa(ii,jj)/24. cm(i,j,13)=6.*gsa(ii,jj)/24. cm(i,j,25)=6.*gsa(ii,jj)/24. cm(i,j,37)=6.*gsa(ii,jj)/24. else cm(i,j,l)=gsa(ii,jj)/24. endif else cm(i,j,l)=0. endif go to 90 endif 90 if(cm(i,j,l).lt.0.) write(99,*) 'Rmorph < 0 at: ', x i,j,ii,jj,l,cmsa(ij,jj),gsa(ii,jj),cm(i,j,l) c if(ida.eq.1) write(99,'(i2,2f7.3,4i5,5f8.3)') l,ylat,xlon, c x i,j,ii,jj,gsa(ii,jj),cmsa(ii,jj),ratio(ii,jj), c x h(l),cm(i,j,l) 100 continue 200 continue 300 continue c return end c c---------------------------------------------------------------- c subroutine adjust_us (cm,icx,jcy,usmask,xlonw,ylatn,gag,cmus, x ratio,ius,jus,thresh,ida) real*4 cm(icx,jcy,48), gag(ius,jus), usmask(ius,jus), h(48), x cmus(ius,jus), ratio(ius,jus) c c c "i1" , 'i2", "j1" & "j2" are indices in the global CMORPH data c ("cm" crray) that correspond to the CPC U.S. rain gauge 1/4 deg. grid c i1=3023 ! 219.875 i2=4122 ! 299.875 j1=1 ! 59.875N j2=621 ! 14.875N Note: CMORPH oriented N-S so enter "-14.875" in Grads to get this value c c c c..... Adjust the half-hourly, 8km CMORPH values by the ratio between c the 1/4 x 1/4 deg gauge & 1/4 x 1/4 deg CMORPH c xbeg=0.036378335 xint=0.072756669 c ybeg=59.963614 yint=0.072771377 do 300 j=j1,j2 ylat=ybeg - (j-1)*yint cc jj=(ylatn-ylat)/0.25 + 1 jj=(ylatn-ylat)/0.25 if(jj.lt.1) jj=1 if(jj.gt.jus) jj=jus do 200 i=i1,i2 xlon=(i-1)*xint + xbeg ii=(xlon-xlonw)/0.25 + 1 if(ii.lt.1) ii=1 if(ii.gt.ius) ii=ius if(usmask(ii,jj).ne.1.) go to 200 if(gag(ii,jj).lt.0.) go to 200 do 100 l=1,48 h(l)=cm(i,j,l) c c.............. If daily CMORPH total is nonzero, then, if the hourly c CMORPH is > zero, multiply each hourly CMORPH value by the ratio c between the daily gauge and daily CMORPH totals c if(cmus(ii,jj).gt.0.) then if(cm(i,j,l).gt.0.) then cm(i,j,l)=cm(i,j,l) * ratio(ii,jj) endif go to 90 endif c c.............. If daily CMORPH total is < or = 0 (i.e. all hourly values are c missing or zero, and the daily gauge total > "thresh", set all c hourly CMORPH values to 1/24th of the daily gauge amt.; otherwise, set c all CMORPH hourly amts. to zero -- this is done because the rain gauge c amounts are small and may be spurious due to the rain c gauge objective analysis technique, i.e. spreading out of low precip c amounts due to smoothing. c if(cmus(ii,jj).le.0.) then if(gag(ii,jj).gt.thresh) then c c................... Because the final RMORPH values are scaled (by "5") c and stored as a scaled integer in a byte, values of c hourly precip that are < 1 will be evaluated as zero. Therefore, if the c daily rain gauge value is < 4.8 (i.e. 24/5), then instead of setting all 24 c hourly values to 1/24th of the daily gauge value for this situation, set c the values of the 4 synoptic times to 1/4th of the daily gauge value. c if(gag(ii,jj).lt.4.8) then cm(i,j,1)= 6.*gag(ii,jj)/24. cm(i,j,13)=6.*gag(ii,jj)/24. cm(i,j,25)=6.*gag(ii,jj)/24. cm(i,j,37)=6.*gag(ii,jj)/24. else cm(i,j,l)=gag(ii,jj)/24. endif else cm(i,j,l)=0. endif go to 90 endif 90 if(cm(i,j,l).lt.0.) write(99,*) 'Rmorph < 0 at: ', x i,j,ii,jj,l,cmus(ij,jj),gag(ii,jj),cm(i,j,l) c if(ida.eq.1) write(99,'(i2,2f7.3,4i5,5f8.3)') l,ylat,xlon, c x i,j,ii,jj,gag(ii,jj),cmus(ii,jj),ratio(ii,jj), c x h(l),cm(i,j,l) 100 continue 200 continue 300 continue c return end c c-------------------------------------------------------------------- c subroutine gauge_sa(io_gauge,iyrmndy,isa,jsa,gsa, x samask,tempsa,iret) real*4 gsa(isa,jsa), samask(isa,jsa), tempsa(isa,jsa) integer*4 jbuf(isa,jsa) character*120 high_level_path, f_gauge logical inthere c c c..... Read in high level directory path from "~cpcsat/.aliases.scripts" file call getarg(2,high_level_path) do i=1,120 if(high_level_path(i:i).eq.' ') go to 110 enddo write(99,*)'WARNING: No blank found in high-level path' 110 len=i-1 c c...... Note: "recent" data (~ past 1 yr) in directory: cwlinks/dev/prec1/data/sa_grid c older data in directory: cwlinks/dev/prec1/data/y----/sa_grid c c f_gauge=high_level_path(1:len) // c x '/cwlinks/dev/prec1/data/sa_grid/y----/yyyymmdd.ll' iyr=iyrmndy/10000 c write (f_gauge(34+len:37+len),'(i4)') iyr c write (f_gauge(39+len:46+len),'(i8)') iyrmndy f_gauge=high_level_path(1:len) // x '/cwlinks/dev/prec1/data/sa_grid/yyyymmdd.ll' write (f_gauge(33+len:40+len),'(i8)') iyrmndy c inquire(file=f_gauge,exist=inthere) if(.not.inthere) then write(99,*) 'NO S. America Gauge file:',f_gauge iret=-1 return endif c c....... Read "data" & orient North > South c write(99,*) 'opening gauge file:',f_gauge open (io_gauge,file=f_gauge,form='unformatted') read(io_gauge) tempsa close(io_gauge) c jj=jsa+1 do j=1,jsa jj=jj-1 do i=1,isa gsa(i,jj)=tempsa(i,j) enddo enddo iret=0 idy=mod(iyrmndy,100) write(99,*) 'exiting gauge_sa' return end c c-------------------------------------------------------------------- c subroutine gauge_us(io_gauge,iyrmndy,ius,jus,gag, x usmask,temp,iret) real*4 gag(ius,jus), usmask(ius,jus), temp(ius,jus) character*120 high_level_path, f_gauge, f_mask logical inthere c c c..... Read in high level directory path from "~cpcsat/.aliases.scripts" file call getarg(3,high_level_path) do i=1,120 if(high_level_path(i:i).eq.' ') go to 110 enddo write(99,*)'WARNING: No blank found in high-level path' 110 len=i-1 c f_gauge=high_level_path(1:len) // '/jj/yyyymmdd.test.ll' write (f_gauge(5+len:12+len),'(i8)') iyrmndy c inquire(file=f_gauge,exist=inthere) if(.not.inthere) then write(99,*) 'NO US Gauge file:',f_gauge iret=-1 return endif c c....... Read "data", convert units from in/day to mm/day & orient North > South c write(99,*) 'opening gauge file:',f_gauge open (io_gauge,file=f_gauge,form='unformatted') read(io_gauge) temp close(io_gauge) c jj=jus+1 do j=1,jus jj=jj-1 do i=1,ius gag(i,jj)=temp(i,j)*25.4 enddo enddo iret=0 write(99,*) 'exiting gauge_us' return end c c---------------------------------------------------------------------- c subroutine grid(iunit,f_gauge,grdint,nx,ny,ntot,gnumbr) real*4 gnumbr(nx,ny) character*120 f_gauge c c------ greates an array that contains the # of rain gauge reports in c each grid box for the S. American daily rain gauge analyses c c open (iunit,file=f_gauge) c xlonw=270. ylatn=15. c c do 200 irec=1,99999 read(iunit,*,err=180,end=250) id,xlon,ylat,prc zlon=xlon if(xlon.lt.0.) zlon=xlon+360. ix = (zlon-xlonw)/grdint + 1 jy = (ylatn-ylat)/grdint + 1 if(ix.gt.nx .or. ix.lt.1 .or. jy.gt.ny .or. jy.lt.1) then write(99,*) xlon,zlon,ylat,ix,jy go to 200 endif gnumbr(ix,jy) = gnumbr(ix,jy) + 1. ntot=ntot+1 go to 200 180 write(99,*) 'Read error at rec # ',irec 200 continue c 250 write(99,*) '# of "raw SA gauge" reports: ',ntot close (iunit) return end c c---------------------------------------------------------------------------- c subroutine gtcmorph(iyr,imo,ida,ihr1,ihr2,icx,jcy,buf,jhr,cm) real*4 cm(icx,jcy,48) character*1 buf(icx,jcy) character*120 fin character*40 jtype character*6 us_technique logical inthere jtype='cmorph_8km_opnl' do 200 ihr=ihr1,ihr2 jhr=jhr+1 call pathfinder(jtype,iyr,imo,ida,ihr,min,idoy,idsat, x us_technique,fin,inthere,iret,len) if(.not.inthere) then write(99,*) 'no CMORPH file:',fin go to 200 endif if(iret.ne.0) stop 'Problem with "pathfinder"' write(99,*) 'opening file:',fin open (2,file=fin,access='direct',recl=icx*jcy*1) c c NOTE: The operational 8km CMORPH files contain 6 records. The 1st c group of 3 records are valid for the 1st half-hour; the 2nd group of c 3 records are valid for the 2nd half hour. THe group of 3 c records are (in order) CMORPH, time-stamp & satellite ID. c do ihalf=1,2 c irec=(ihalf-1)*3+1 read(2,rec=irec) buf c khr=(jhr-1)*2+ihalf write(99,*) 'jhr, khr: ',jhr,khr do j=1,jcy do i=1,icx ival=ichar(buf(i,j)) if(ival.eq.255) then cm(i,j,khr)=-9999. else c................. if .not. missing, then unscale to get mm/hr units cm(i,j,khr)=0.2*ival endif enddo enddo enddo close(2) 200 continue return end c c---------------------------------------------------------------------- c subroutine gtname(iyr,imo,ida,f_gauge,iret) c c----- Builds the pathname of the S. American daily raingauge analysis c for the specified date c c character*120 high_level_path, f_gauge logical inthere c c..... Read in high level directory path from "~cpcsat/.aliases.scripts" file call getarg(2,high_level_path) do i=1,120 if(high_level_path(i:i).eq.' ') go to 110 enddo write(99,*)'WARNING: No blank found in high-level path' 110 len=i-1 c f_gauge=high_level_path(1:len) // x'/cwlinks/dev/prec1/data/sa_raw/raw_data/y----/prec0d0myyyy12.dat' c/cwlinks/dev/prec1/data/sa_raw/raw_data/y----/precddmmyyyy12.dat write (f_gauge(42+len:45+len),'(i4)') iyr write (f_gauge(51+len:52+len),'(i2)') ida if(ida.lt.10) f_gauge(51+len:51+len)='0' write (f_gauge(53+len:54+len),'(i2)') imo if(imo.lt.10) f_gauge(53+len:53+len)='0' write (f_gauge(55+len:58+len),'(i4)') iyr c inquire(file=f_gauge,exist=inthere) if(.not.inthere) then write(99,*) 'NO Gauge file:',f_gauge iret=-1 else write(99,*) 'opening gauge "raw data" file:',f_gauge iret=1 endif return end c c-------------------------------------------------------------------- c c subroutine masksa(io_mask,isa,jsa,tempsa,samask) character*120 high_level_path, f_mask real*4 samask(isa,jsa), tempsa(isa,jsa), omask(81,101) c c..... Read in high level directory path from "~cpcsat/.aliases.scripts" file call getarg(2,high_level_path) do i=1,120 if(high_level_path(i:i).eq.' ') go to 110 enddo write(99,*)'WARNING: No blank found in high-level path' 110 len=i-1 c f_mask=high_level_path(1:len) // x '/cwlinks/dev/prec1/data/sa_grid/sa_mask1x1.dat' write(99,*) 'opening mask file:',f_mask open (io_mask,file=f_mask,form='unformatted') read(io_mask) omask close(io_mask) c c......... SA mask is oriented North->South so no need to c reverse array jj=0 do j=6,81 jj=jj+1 ii=0 do i=21,81 ii=ii+1 tempsa(ii,jj)=omask(i,j) enddo enddo do j=1,jsa do i=1,isa samask(i,j)=tempsa(i,j) enddo enddo c...... Because mask 'bleeds' into surrounding ocean on northern c part of SA continent, modify some mask values to c indicate that grids are ocean covered rather than land c samask(29,1)=-1. samask(9,2)=-1. samask(10,2)=-1. do i=29,32 samask(i,2)=-1. samask(i,3)=-1. enddo samask(8,3)=-1. samask(9,3)=-1. do i=19,23 samask(i,3)=-1. enddo samask(8,4)=-1. do i=16,18 samask(i,4)=-1. enddo do i=23,30 samask(i,4)=-1. samask(i,5)=-1. enddo samask(4,5)=-1. samask(8,5)=-1. do i=9,15 samask(i,6)=-1. enddo samask(6,7)=-1. samask(31,7)=-1. samask(32,7)=-1. samask(7,8)=-1. samask(12,8)=-1. samask(34,9)=-1. return end c c-------------------------------------------------------------------- c c subroutine maskus(io_mask,ius,jus,usmask) character*120 high_level_path, f_mask real*4 usmask(ius,jus) c c..... Read in high level directory path from "~cpcsat/.aliases.scripts" file call getarg(2,high_level_path) do i=1,120 if(high_level_path(i:i).eq.' ') go to 110 enddo write(99,*)'WARNING: No blank found in high-level path' 110 len=i-1 c f_mask=high_level_path(1:len) // x '/cwlinks/dev/prec1/data/us_mex_grid/us_mex_mask0.25x0.25.dat' write(99,*) 'opening mask file:',f_mask open (io_mask,file=f_mask,form='unformatted') read(io_mask) usmask close(io_mask) c c..... US mask is oriented N-> S so do not reverse c return end c c-------------------------------------------------------------------- c subroutine mk_ratio_sa(ylatn,xlonw,cm,icx,jcy,samask,gsa,cmsa, x xpts,ratio,isa,jsa,nday) c c................................................................................ c Purpose: Accumulates half-hourly CMORPH totals to daily totals c that matches the daily accumulation period for S. AMerica (12Z-12Z) c for 1 x 1 degree grids that match the CPC S. American gauge analysis. c Then routine computes the ratio between the daily c CMORPH & daily gauge data which is the applied (in subroutine c "adjust_sa") to make the RMORPH values. c............................................................................... c real*4 cm(icx,jcy,48), gsa(isa,jsa), cmsa(isa,jsa), x samask(isa,jsa), xpts(isa,jsa), ratio(isa,jsa), h(48) c do j=1,jsa do i=1,isa cmsa(i,j)=0. xpts(i,j)=0. enddo enddo c c c "ii" & "jj" are indices that map the Cmorph data to the c appropriate S. American analysis gridbox. c i1=3713 ! 270.125 i2=4538 ! 330.125 c j1= 621 ! 14.875N Note: CMORPH oriented N-S so enter "-14.875" in Grads to get this value j2=1649 ! 59.875S c c xbeg=0.036378335 xint=0.072756669 c ybeg=59.963614 yint=0.072771377 c c......Compute 24-h 1x1 deg. CMORPH areal mean (to match gauge spatial resolution) c do j=j1,j2 ylat=ybeg - (j-1)*yint jj=ylatn-ylat+1. write(99,*) j,ylat,jj do i=i1,i2 xlon=xbeg+xint*(i-1) ii=xlon-xlonw + 1 if(j.eq.j1) write(99,*) i,xlon,ii ccc write(99,*) ylat,xlon,jj,ii xpts(ii,jj)=xpts(ii,jj)+1. do l=1,48 if(cm(i,j,l).ge.0.) cmsa(ii,jj)=cmsa(ii,jj) + cm(i,j,l) enddo enddo enddo c do j=1,jsa do i=1,isa if(xpts(i,j).gt.0.)then cmsa(i,j)=cmsa(i,j)/xpts(i,j) else cmsa(i,j)=-9999. endif c c...... Compute the ratio between the gauge analysis & CMORPH (at the gauge spatial scale) c if(cmsa(i,j).ge.0.5) then ratio(i,j)=gsa(i,j)/cmsa(i,j) else ratio(i,j)=1.0 endif if(samask(i,j).ne.1.) ratio(i,j)=-9999. c write(99,*) i,j,gsa(i,j), cmsa(i,j), ratio(i,j) enddo enddo write(98,rec=(nday-1)*6+3) cmsa write(98,rec=(nday-1)*6+4) ratio write(98,rec=(nday-1)*6+5) xpts return end c c-------------------------------------------------------------------- c subroutine mk_ratio_us(ylatn,xlonw,cm,icx,jcy,usmask,gag,cmus, x xpts,ratio,ius,jus,nday) c c................................................................................ c Purpose: Accumulates half-hourly CMORPH totals to daily totals c that matches the daily accumulation period for N. America (12Z-12Z) c for 1/4 x 1/4 degree grids that match the CPC N. American gauge analysis. c Then routine computes the ratio between the daily c CMORPH & daily gauge data which is the applied (in subroutine c "adjust_us") to make the RMORPH values. c............................................................................... c real*4 cm(icx,jcy,48), gag(ius,jus), cmus(ius,jus), x usmask(ius,jus), xpts(ius,jus), ratio(ius,jus), h(48) c do j=1,jus do i=1,ius cmus(i,j)=0. xpts(i,j)=0. enddo enddo c c c "ii" & "jj" are indices that map the Cmorph data to the c appropriate US raingauge analysis gridbox. c i1=3023 ! 219.875 i2=4122 ! 299.875 c j1= 1 ! 59.875N j2= 621 ! 14.875N Note: CMORPH oriented N-S so enter "-14.875" in Grads to get this value c c xbeg=0.036378335 xint=0.072756669 c ybeg=59.963614 yint=0.072771377 c c......Compute 24-h 1/4 x 1/4 deg. CMORPH areal mean (to match gauge spatial resolution) c do j=j1,j2 ylat=ybeg - (j-1)*yint cc jj=(ylatn-ylat)/0.25 + 1 jj=(ylatn-ylat)/0.25 if(jj.lt.1) jj=1 if(jj.gt.jus) jj=jus do i=i1,i2 xlon=xbeg+xint*(i-1) ii=(xlon-xlonw)/0.25 + 1 if(ii.lt.1) ii=1 if(ii.gt.ius) ii=ius ccc write(99,*) '"mk_ratio_us:',xlon,ylat,ii,jj xpts(ii,jj)=xpts(ii,jj)+1. do l=1,48 if(cm(i,j,l).ge.0.) cmus(ii,jj)=cmus(ii,jj) + cm(i,j,l) enddo enddo enddo c do j=1,jus do i=1,ius if(xpts(i,j).gt.0.)then cmus(i,j)=cmus(i,j)/xpts(i,j) else cmus(i,j)=-9999. endif c c...... Compute the ratio between the gauge analysis & CMORPH (at the gauge spatial scale) c if(cmus(i,j).ge.0.5) then ratio(i,j)=gag(i,j)/cmus(i,j) else ratio(i,j)=1.0 endif if(usmask(i,j).ne.1.) ratio(i,j)=-9999. enddo enddo write(78,rec=(nday-1)*4+2) cmus write(78,rec=(nday-1)*4+3) ratio write(78,rec=(nday-1)*4+4) xpts return end c c-------------------------------------------------------------------- c subroutine smth9(datin,gagnum,nx,ny,datout) c c Subroutine does 9-pt smoothing on input data. Each of the 9 c grids is weighted by the number of rain gauge observations c in that grid location. c c Inputs: c c datin: input data to be smoothed c gagnum: array that contains # gauge obs in each grid box c nx, ny: array dimensions c c Oututs: c c datout: array containing smoothed data c c c Programmer: Janowiak April 8, 2006 c c--------------------------------------------------------------------------------- c real*4 datin(nx,ny), datout(nx,ny), gagnum(nx,ny) do j=1,ny do i=1,nx sum=0. wgts=0. jfrst=j-1 jlast=j+1 if(jfrst.le.0) jfrst=1 if(jlast.gt.ny) jlast=ny do jj=jfrst,jlast ifrst=i-1 ilast=i+1 if(ifrst.le.0) ifrst=1 if(ilast.gt.nx) ilast=nx do ii=ifrst,ilast if(datin(ii,jj).ge.0.) then sum = sum + datin(ii,jj)*gagnum(ii,jj) wgts = wgts + gagnum(ii,jj) endif enddo enddo if(wgts.gt.0.) then datout(i,j) = sum/wgts else datout(i,j)=datin(i,j) endif enddo enddo return end c c-------------------------------------------------------------------- c subroutine wtrmorph(iyrlast,imolast,idalast,iyr,imo,ida,icx,jcy, x buf,cm) real*4 cm(icx,jcy,48) character*1 buf(icx,jcy) character*120 outfl data imsng /255/ outfl='data/yyyy/yyyymmddhh' c c c Writes out the "RMORPH" data. c c NOTE: subroutine is set up for 12z-12z c periods, so modifications will be necessary for daily rain gauge data sets that c have different "day definitions" c c c Data in "cm" array are stored in the following order: c c 1: 12:00z c 2: 12:30z c ... c ..........(Next day) c 25: 00:00z c 26: 00:30z c ... c 47: 11:00z c 48: 11:30z c c So, "jhr" is the correct hour for naming output files c do 200 ihr=1,12 jhr=ihr+11 iyrmndyhr=iyrlast*1000000+imolast*10000+idalast*100+jhr write(outfl(6:9),'(i4)') iyrlast write(outfl(11:20),'(i10)') iyrmndyhr write(99,*) 'opening output file:',outfl open(10,file=outfl,access='direct',recl=icx*jcy*1) do ihalf=1,2 khr=(ihr-1)*2 + ihalf write(99,*) 'wtrmorph (ihr,jhr,khr)', ihr,jhr,khr do j=1,jcy do i=1,icx if(cm(i,j,khr).eq.-9999.) then buf(i,j)=char(imsng) else c................. if .not. missing, then scale ival= 5. * cm(i,j,khr) buf(i,j)=char(ival) c if (i.eq.956 .and. j.eq.375) write(99,*) 'DISAG:',ihr,ida, c x cm(i,j,ihr),ival endif enddo enddo write(10,rec=ihalf) buf enddo close(10) 200 continue c do 300 ihr=13,24 jhr=ihr-13 iyrmndyhr=iyr*1000000+imo*10000+ida*100+jhr write(outfl(6:9),'(i4)') iyr write(outfl(11:20),'(i10)') iyrmndyhr write(99,*) 'opening output file:',outfl open(10,file=outfl,access='direct',recl=icx*jcy*1) do ihalf=1,2 khr=(ihr-1)*2 + ihalf write(99,*) 'wtrmorph (ihr,jhr,khr)', ihr,jhr,khr do j=1,jcy do i=1,icx if(cm(i,j,khr).eq.-9999.) then buf(i,j)=char(imsng) else c................. if .not. missing, then scale ival= 5. * cm(i,j,khr) buf(i,j)=char(ival) c if (i.eq.956 .and. j.eq.375) write(99,*) 'DISAG:',ihr,ida, c x cm(i,j,ihr),ival endif enddo enddo write(10,rec=ihalf) buf enddo close(10) 300 continue 500 continue return end