subroutine spruv(wdatain,nn_w,ugesin,vgesin, & 1,3 fact,wtype,nwdta,nsig,ermax,ermin,gross,lat12,& nx,nlon,weight,mype,awork,wpres,nlev,ptop,pbot,vwork) !$$$ subprogram documentation block ! . . . . ! subprogram: spruv save wind infor. ! prgmmr: parrish org: w/nmc22 date: 90-10-13 ! ! abstract: save wind information for later use. ! ! program history log: ! 90-10-13 parrish ! 98-04-07 weiyu yang ! 99-08-24 derber, j., treadon, r., yang, w., first frozen mpp version ! ! input argument list: ! wdatain - obs info at obs locations in each pe sub-domain ! nn_w - number of data in each wind record ! fact - reduction to 10,20 m factor ! wtype - wind type ! nwdta - number of obs ! nsig - number of layers on gaussian grid ! ermax,ermin,gross - parameters for gross error test ! lat12 - lat1+2 ! nlon - number of longitudes ! weight - location weighting ! mype - pe number ! ! output argument list: ! j_w - wind data location indices ! wgt_w - weighting information of wind data ! other_w - other information of wind data, other_w ! is multiplied by aerr already. ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use type_kinds, only: fp_kind use obsmod implicit none include 'mpi_inc.h' include 'ijstart.h' integer iout_uv,mype_uv integer jlat,jlon,jsig,ict,icts,m1,n1,n2,ier,ilon,ilat,isig integer iures,ivres,iii,k,ioff,i,lat12,nn_w,nwdta,nsig,iix integer itypx,itype,latlon11_l,nlev,nx,nlon,mype real(fp_kind) dx,dy,ds,dx1,dy1,ds1,scale,ratio,obserror,obserrlm,small real(fp_kind) residual,val,ressw,ress,val1,val2,ermin,gross,ermax real(fp_kind),dimension(nwdta,nn_w):: wdatain real(fp_kind),dimension(nwdta):: wtype,fact,ugesin,vgesin real(fp_kind),dimension(nwdta):: weight,wpres real(fp_kind),dimension(100+7*nsig,6):: awork real(fp_kind),dimension(6):: ptop,pbot real(fp_kind),dimension(6,101,3):: vwork small=1.e-10 nwdata_dim=max(nwdata_dim,1) nspddat_dim=max(nspddat_dim,1) call makew m1=mype+1; n1=nsig+3; n2=n1+nsig ier=1; ilon=2; ilat=3; isig=4; iures=5; ivres=6 do iii=1,nwdta if(wdatain(iii,isig)<1.0) wdatain(iii,isig)=1.0 end do ict = 0 icts = 0 do iix=1,2 do iii=1,nwdta itypx=nint(wtype(iii)) if(itypx/=283.and.iix==1)then jlat=wdatain(iii,ilat) jlon=wdatain(iii,ilon) jsig=wdatain(iii,isig) dx=wdatain(iii,ilon)-jlon; dx1=1.0-dx dy=wdatain(iii,ilat)-jlat; dy1=1.0-dy ds=wdatain(iii,isig)-jsig; ds1=1.0-ds jlat=min(max(1,jlat),nx); jlon=min(max(0,jlon),nlon) latlon11_l=latlon11_obs if(jsig==nsig) latlon11_l=0 if(jstart(m1)==1.and.jlon==nlon) jlon=0 if(jstart(m1)+jlon1(m1)==nlon+1.and.jlon==0) jlon=nlon jlat=jlat-istart(m1)+2; jlon=jlon-jstart(m1)+2 !-----------------------------------gross error test added here obserror=1.0/max(sqrt(wdatain(iii,ier)),small) obserrlm=max(ermin,min(ermax,obserror)) residual=sqrt(wdatain(iii,iures)**2+wdatain(iii,ivres)**2) ratio=residual/obserrlm if(ratio>gross) then awork(4,1)=awork(4,1)+weight(iii); wdatain(iii,ier)=0.0 end if awork(4*nsig+jsig+100,1)=awork(4*nsig+jsig+100,1)+wdatain(iii,ier)*& wdatain(iii,iures)**2*weight(iii) awork(5*nsig+jsig+100,1)=awork(5*nsig+jsig+100,1)+wdatain(iii,ier)*& wdatain(iii,ivres)**2*weight(iii) awork(6*nsig+jsig+100,1)=awork(6*nsig+jsig+100,1)+weight(iii) wdatain(iii,ier)=sqrt(wdatain(iii,ier)) if(wdatain(iii,ier) > 1.e-10)then ict=ict+1 j1_w(ict)=jlat+(jlon-1)*lat12+(jsig-1)*latlon11_obs j2_w(ict)=j1_w(ict)+1 j3_w(ict)=j1_w(ict)+lat12 j4_w(ict)=j3_w(ict)+1 j5_w(ict)=j1_w(ict)+latlon11_l j6_w(ict)=j5_w(ict)+1 j7_w(ict)=j5_w(ict)+lat12 j8_w(ict)=j7_w(ict)+1 w1_w(ict)=fact(iii)*wdatain(iii,ier)*dx1*dy1*ds1 !wgt000 w2_w(ict)=fact(iii)*wdatain(iii,ier)*dx1*dy*ds1 !wgt100 w3_w(ict)=fact(iii)*wdatain(iii,ier)*dx*dy1*ds1 !wgt010 w4_w(ict)=fact(iii)*wdatain(iii,ier)*dx*dy*ds1 !wgt110 w5_w(ict)=fact(iii)*wdatain(iii,ier)*dx1*dy1*ds !wgt001 w6_w(ict)=fact(iii)*wdatain(iii,ier)*dx1*dy*ds !wgt101 w7_w(ict)=fact(iii)*wdatain(iii,ier)*dx*dy1*ds !wgt011 w8_w(ict)=fact(iii)*wdatain(iii,ier)*dx*dy*ds !wgt111 ures(ict)=wdatain(iii,iures)*wdatain(iii,ier) vres(ict)=wdatain(iii,ivres)*wdatain(iii,ier) uges(ict)=ugesin(iii)*wdatain(iii,ier) vges(ict)=vgesin(iii)*wdatain(iii,ier) uverr(ict)=wdatain(iii,ier) weight_w(ict)=weight(iii) end if end if if(itypx==283.and.iix==2)then jlat=wdatain(iii,ilat) jlon=wdatain(iii,ilon) dx=wdatain(iii,ilon)-jlon; dx1=1.0-dx dy=wdatain(iii,ilat)-jlat; dy1=1.0-dy if(jstart(m1)==1.and.jlon==nlon) jlon=0 if(jstart(m1)+jlon1(m1)==nlon+1.and.jlon==0) jlon=nlon jlat=jlat-istart(m1)+2; jlon=jlon-jstart(m1)+2 !-----------------------------------gross error test added here obserror=1.0/max(sqrt(wdatain(iii,ier)),small) obserrlm=max(ermin,min(ermax,obserror)) residual=abs(wdatain(iii,iures)-wdatain(iii,ivres)) ratio=residual/obserrlm if(ratio>gross) then awork(4,1)=awork(4,1)+weight(iii); wdatain(iii,ier)=0.0 end if awork(5,1)=awork(5,1)+wdatain(iii,ier)*& (wdatain(iii,iures)-wdatain(iii,ivres))**2*weight(iii) awork(6,1)=awork(6,1)+weight(iii) wdatain(iii,ier)=sqrt(wdatain(iii,ier)) if(wdatain(iii,ier) > 1.e-10)then icts=icts+1 j1_spd(icts)=jlat+(jlon-1)*lat12 j2_spd(icts)=j1_spd(icts)+1 j3_spd(icts)=j1_spd(icts)+lat12 j4_spd(icts)=j3_spd(icts)+1 w1_spd(icts)=fact(iii)*wdatain(iii,ier)*dx1*dy1 !wgt000 w2_spd(icts)=fact(iii)*wdatain(iii,ier)*dx1*dy !wgt100 w3_spd(icts)=fact(iii)*wdatain(iii,ier)*dx*dy1 !wgt010 w4_spd(icts)=fact(iii)*wdatain(iii,ier)*dx*dy !wgt110 spdres(icts)=wdatain(iii,iures)*wdatain(iii,ier) usges(icts)=ugesin(iii)*wdatain(iii,ier) vsges(icts)=vgesin(iii)*wdatain(iii,ier) spderr(icts)=wdatain(iii,ier) weight_spd(icts)=weight(iii) end if end if end do end do nwdata_num=ict if(nwdata_num /= nwdata_dim)print *,' wind dim, num ',nwdata_dim,nwdata_num nspddat_num=icts if(nspddat_num /= nspddat_dim)print *,' spd dim, num ',nspddat_dim,nspddat_num ioff=200 scale=1.0 do k = 1,nlev do i=1,nwdta if (wpres(i)>=ptop(k).and.wpres(i)<=pbot(k)) then val1=wdatain(i,ier)*wdatain(i,iures) val2=wdatain(i,ier)*wdatain(i,ivres) if (nint(wdatain(i,9))/=283) then ress=scale*sqrt(wdatain(i,5)**2+wdatain(i,6)**2) val=.5*(val1*val1+val2*val2) else ress=scale*(wdatain(i,5)-wdatain(i,6)) val=(val1-val2)**2 end if val=weight(i)*val ressw=ress*ress*weight(i) itype=nint(wdatain(i,9))-ioff itype=max(1,min(itype,100)) vwork(k,itype,1) = vwork(k,itype,1)+weight(i) vwork(k,itype,2) = vwork(k,itype,2)+ressw vwork(k,itype,3) = vwork(k,itype,3)+val vwork(k,101,1) = vwork(k,101,1)+weight(i) vwork(k,101,2) = vwork(k,101,2)+ressw vwork(k,101,3) = vwork(k,101,3)+val end if end do end do return end