cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c subroutine LLWS: compute low level wind shear c c Author: Binbin Zhou, Aug, 7, 2005 c 05-17-2013, B. Zhou: Upgraded to grib2 I/O c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine llws (nv,ifunit,jf,iens,Lm,Lp,Lt,eps, + derv_mn,derv_sp,derv_pr,wgt) use grib_mod include 'parm.inc' parameter (lvl=10) !max 10 levels is set to reduce memory !parameter (lvl=14) !max 14 levels is set c for derived variables Character*4 dvname(maxvar) Integer dk5(maxvar), dk6(maxvar), dk4(maxvar) Character*1 dMsignal(maxvar), dPsignal(maxvar) Integer dMlvl(maxvar), dMeanLevel(maxvar,maxmlvl) Integer dPlvl(maxvar), dProbLevel(maxvar,maxplvl) Character*1 dop(maxvar) Integer dTlvl(maxvar) Real dThrs(maxvar,maxtlvl) Integer MPairLevel(maxvar,maxmlvl,2) Integer PPairLevel(maxvar,maxplvl,2) common /dtbl/nderiv, + dvname,dk4,dk5,dk6,dMlvl,dPlvl,dTlvl, + dMeanLevel,dProbLevel,dThrs, + dMsignal,dPsignal,MPairLevel,PPairLevel,dop INTEGER, intent(IN) :: nv, jf, iens REAL,dimension(jf,Lm),intent(INOUT) :: derv_mn,derv_sp REAL,dimension(jf,Lm,Lt),intent(INOUT) :: derv_pr real apoint(iens),LLWSapoint(iens),llws1(jf,iens) real ws,u10,v10,hsfc,up(lvl),vp(lvl),hp(lvl) INTEGER miss(iens) real wgt(30) integer,dimension(iens),intent(IN) :: ifunit type(gribfield) :: gfld REAL, dimension(jf,iens)::U10m,V10m,HS REAL, dimension(jf,iens,lvl)::Upr,Vpr,Hpr Integer p(50),lvls character*5 eps data (p(i),i=1,25) + /1000,975,950,925,900,875,850,825,800,775,750,725,700,675, + 650,625,600,575,550,525,500,475,450,425,400/ !write(*,*) 'In LLWS .....' !write(*,*) 'ifunit=', ifunit, eps lvls=10 write(*,*) nv,jf,iens,Lm,Lp,Lt,lvls jpdtn=0 jp27=-9999 miss=0 do irun=1,iens call readGB2(ifunit(irun),jpdtn,2,2,103,10,jp27,gfld,ie) !U10m U10m(:,irun)=gfld%fld call readGB2(ifunit(irun),jpdtn,2,3,103,10,jp27,gfld,ie) !V10m V10m(:,irun)=gfld%fld call readGB2(ifunit(irun),jpdtn,3,5,1,0,jp27,gfld,ie) !Hsfc HS(:,irun)=gfld%fld !write(*,*) 'read surface data done' do k=1,lvls !write(*,*) 'Level=',p(k) call readGB2(ifunit(irun),jpdtn,2,2,100,p(k),jp27,gfld,ie) Upr(:,irun,k)=gfld%fld call readGB2(ifunit(irun),jpdtn,2,3,100,p(k),jp27,gfld,ie) Vpr(:,irun,k)=gfld%fld call readGB2(ifunit(irun),jpdtn,3,5,100,p(k),jp27,gfld,ie) Hpr(:,irun,k)=gfld%fld end do end do do lv=1,dMlvl(nv) !for all levels do igrid = 1,jf do i=1,iens u10=U10m(igrid,i) v10=V10m(igrid,i) hsfc=HS(igrid,i) do k=1,lvls up(k)=Upr(igrid,i,k) vp(k)=Vpr(igrid,i,k) hp(k)=Hpr(igrid,i,k) end do if (u10.eq.0.0.and.v10.eq.0.0.and.hsfc.eq.0.0) then ws=0.0 goto 500 end if call get_llws(up,vp,hp,u10,v10,hsfc,jf,lvls,ws) 500 LLWSapoint(i)=ws llws1(igrid,i)=ws end do call getmean(LLWSapoint,iens,amean,aspread, + miss,wgt) derv_mn(igrid,lv)=amean derv_sp(igrid,lv)=aspread c if(amean.gt.20.0) write(*,*) 'LLWS:',igrid,amean !if(igrid.eq.10000) then ! write(*,*) 'mean LLWS=',amean,aspread !end if end do end do do lv=1,dPlvl(nv) do lt = 1, dTlvl(nv) do igrid = 1,jf LLWSapoint(:)=llws1(igrid,:) if(trim(dop(nv)).ne.'-') then thr1 = dThrs(nv,lt) thr2 = 0. call getprob(LLWSapoint,iens,thr1,thr2,dop(nv),aprob, + miss,wgt) derv_pr(igrid,lv,lt)=aprob else if(lt.lt.dTlvl(nv)) then thr1 = dThrs(nv,lt) thr2 = dThrs(nv,lt+1) call getprob(LLWSapoint,iens,thr1,thr2,dop(nv),aprob, + miss,wgt) derv_pr(igrid,lv,lt)=aprob end if end if end do end do end do return end subroutine get_llws (up,vp,hp,u10,v10,hsfc,jf,lv,llws) c up,vp : u,v on pressure levels c hp : height of pressure levels c u10,v10 : 10m u and v c hsfc : sfc height INTEGER jf,lv REAL up(lv),vp(lv),hp(lv) REAL u10,v10,hsfc,llws !write(*,*) 'in get_llws ...' !write(*,*) u10,v10,hsfc,jf,lv !do k=1,lv ! write(*,*)k, up(k),vp(k),hp(k) !end do z1 = 10. + hsfc if (z1 .lt. hp(1) ) then kp = 0 else do k = 1,lv-1 if(z1.ge.hp(k).and.z1.lt.hp(k+1)) then kp=k end if end do end if hz1 = hp(kp+1) - z1 dh=0.0 if((hz1+10.0).gt.609.6) then !609.6m=2000ft, hz1+10 means from sfc instead form 10m u2=u10 + (up(kp+1)-u10)*599.6/hz1 v2=v10 + (vp(kp+1)-v10)*599.6/hz1 z2=hsfc+609.6 else do k=kp+1,lv-1 dh=dh+hp(k+1)-hp(k) if((dh+hz1+10).gt.609.6) then !ie, reach to 2000 feet z2=hsfc+609.6 rt=(z2-hp(k))/(hp(k+1)-hp(k)) u2=up(k)+(up(k+1)-up(k))*rt v2=vp(k)+(vp(k+1)-vp(k))*rt k2=k goto 609 end if end do end if 609 llws=abs(sqrt((u2-u10)**2+(v2-v10)**2)) & /609.9 llws=llws*1.943*609.9 !unit--> knot/2000ft return end