! SUBSET=NC006001 -- level 3 superobs ! SUBSET=NC006002 -- level 2.5 superobs ! SUBSET=NC006070 -- RADIAL WIND FROM P3 RADAR ! SUBSET=NC0062XX -- RADIAL WIND FROM Level 2 radar obs on grid-tilt(column-tilt) subroutine read_radar(nread,ndata,nodata,infile,lunout,obstype,twind,sis,hgtl_full,nobs) !$$$ subprogram documentation block ! . . . . ! subprogram: read_radar read radar radial winds ! prgmmr: yang org: np23 date: 1998-05-15 ! ! abstract: This routine reads radar radial wind files. ! ! When running the gsi in regional mode, the code only ! retains those observations that fall within the regional ! domain ! ! program history log: ! 1998-05-15 yang, weiyu ! 1999-08-24 derber, j., treadon, r., yang, w., first frozen mpp version ! 2004-06-16 treadon - update documentation ! 2004-07-29 treadon - add only to module use, add intent in/out ! 2005-06-10 devenyi/treadon - correct subset declaration ! 2005-08-02 derber - modify to use convinfo file ! 2005-09-08 derber - modify to use input group time window ! 2005-10-11 treadon - change convinfo read to free format ! 2005-10-17 treadon - add grid and earth relative obs location to output file ! 2005-10-18 treadon - remove array obs_load and call to sumload ! 2005-10-26 treadon - add routine tag to convinfo printout ! 2006-02-03 derber - modify for new obs control and obs count ! 2006-02-08 derber - modify to use new convinfo module ! 2006-02-24 derber - modify to take advantage of convinfo module ! 2006-04-21 parrish - modify to use level 2, 2.5, and/or 3 radar wind ! superobs, with qc based on vad wind data. ! 2006-05-23 parrish - interpolate model elevation to vad wind site ! 2006-07-28 derber - use r1000 from constants ! 2007-03-01 tremolet - measure time from beginning of assimilation window ! 2008-04-17 safford - rm unused vars and uses ! 2008-09-08 lueken - merged ed's changes into q1fy09 code ! 2009-06-08 parrish - remove erroneous call to cosd, sind ! 2009-05-08 tong - add reading NOAA P3 tail Dopple radar data ! 2010-09-08 parrish - remove subroutine check_rotate_wind. This was a debug routine introduced when ! the reference wind rotation angle was stored as an angle, beta_ref. This field ! had a discontinuity at the date line (180E), which resulted in erroneous wind ! rotation angles for a small number of winds whose rotation angle was interpolated ! from beta_ref values across the discontinuity. This was fixed by replacing the ! beta_ref field with cos_beta_ref, sin_beta_ref. ! 2011-03-28 s.liu - add subtype to radial wind observation and limit the use ! of level2.5 and level3 data in Conus domain for NMM and NMMB ! 2011-08-01 lueken - remove deter_zsfc_model (placed in deter_sfc_mod) and fix indentation ! 2012-01-11 m.Hu - add subtype to radial wind observation and limit the use ! of level2.5 and level3 data in Conus domain for ARW ! 2012-06-26 y.li/x.wang add TDR fore/aft sweep separation for thinning,xuguang.wang@ou.edu ! 2012-04-28 s.liu - use new VAD wind ! 2012-11-12 s.liu - add new VAD wind flag ! 2013-01-26 parrish - change from grdcrd to grdcrd1 (to allow successful debug compile on WCOSS) ! 2013-05-07 tong - add reading tdr superobs data ! 2013-05-22 tong - Modified the criteria of seperating fore and aft sweeps for TDR NOAA/FRENCH antenna ! 2015-02-23 Rancic/Thomas - add thin4d to time window logical ! 2015-10-01 guo - consolidate use of ob location (in deg) ! 2016-09-xx CAPS(G. Zhao) - add capability reading CAPS column-tilt radar radial wind data ! 2016-12-21 lippi/carley - add logic to run l2rw loop (==0) or run loop for l3rw and l2_5rw (==1,2) ! to help fix a multiple data read bug (when l2rwbufr and radarbufr were both ! listed in the OBS_INPUT table) and for added flexibility for experimental setups. ! 2018-02-15 wu - add code for fv3_regional option ! 2020-05-04 wu - no rotate_wind for fv3_regional ! ! ! input argument list: ! infile - file from which to read BUFR data ! lunout - unit to which to write data for further processing ! obstype - observation type to process ! twind - input group time window (hours) ! hgtl_full- 3d geopotential height on full domain grid ! ! output argument list: ! nread - number of doppler lidar wind observations read ! ndata - number of doppler lidar wind profiles retained for further processing ! nodata - number of doppler lidar wind observations retained for further processing ! sis - satellite/instrument/sensor indicator ! nobs - array of observations on each subdomain for each processor ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block use kinds, only: r_kind,r_single,r_double,i_kind,i_byte use constants, only: zero,zero_single,half,one,two,three,deg2rad,rearth,rad2deg, & one_tenth,r10,r1000,r60inv,r100,r400,grav_equator, & eccentricity,somigliana,grav_ratio,grav, & semi_major_axis,flattening,two use qcmod, only: erradar_inflate,vadfile,newvad use obsmod, only: iadate,ianldate,l_foreaft_thin use gsi_4dvar, only: l4dvar,l4densvar,iwinbgn,winlen,time_4dvar,thin4d use gridmod, only: regional,nlat,nlon,tll2xy,rlats,rlons,rotate_wind_ll2xy,nsig,& fv3_regional use gridmod, only: wrf_nmm_regional,nems_nmmb_regional,cmaq_regional,wrf_mass_regional use gridmod, only: fv3_regional use convinfo, only: nconvtype,ctwind, & ncmiter,ncgroup,ncnumgrp,icuse,ictype,ioctype,ithin_conv,rmesh_conv,pmesh_conv use convthin, only: make3grids,map3grids,del3grids,use_all use deter_sfc_mod, only: deter_sfc2,deter_zsfc_model use mpimod, only: npe use gsi_io, only: verbose use mpimod, only: mype use directDA_radaruse_mod, only: oe_rw, lvldbg, refl_lowbnd_rw use directDA_radaruse_mod, only: l_correct_azmu, l_correct_tilt, i_correct_tilt, & l_azm_east1st, l_plt_diag_rw use directDA_radaruse_mod, only: l_use_rw_columntilt implicit none ! Declare passed variables character(len=*),intent(in ) :: obstype,infile character(len=20),intent(in ) :: sis real(r_kind) ,intent(in ) :: twind integer(i_kind) ,intent(in ) :: lunout integer(i_kind) ,intent(inout) :: nread,ndata,nodata integer(i_kind),dimension(npe) ,intent(inout) :: nobs real(r_kind),dimension(nlat,nlon,nsig),intent(in):: hgtl_full ! Declare local parameters integer(i_kind),parameter:: maxlevs=1500 integer(i_kind),parameter:: maxdat=22 integer(i_kind),parameter:: maxvad=500 ! integer(i_kind),parameter:: maxvadbins=20 integer(i_kind),parameter:: maxvadbins=15 real(r_kind),parameter:: r4_r_kind = 4.0_r_kind real(r_kind),parameter:: dzvad=304.8_r_kind ! vad reports are every 1000 ft = 304.8 meters real(r_kind),parameter:: r3_5 = 3.5_r_kind real(r_kind),parameter:: r6 = 6.0_r_kind real(r_kind),parameter:: r8 = 8.0_r_kind real(r_kind),parameter:: r90 = 90.0_r_kind real(r_kind),parameter:: r200 = 200.0_r_kind real(r_kind),parameter:: r150 = 150.0_r_kind real(r_kind),parameter:: r360=360.0_r_kind real(r_kind),parameter:: r50000 = 50000.0_r_kind real(r_kind),parameter:: r60 = 60.0_r_kind real(r_kind),parameter:: r75 = 75.0_r_kind real(r_kind),parameter:: r92 = 92.6e03_r_kind real(r_kind),parameter:: r89_5 = 89.5_r_kind real(r_kind),parameter:: r2 = 2.0_r_kind real(r_kind),parameter:: r71 = 71.0_r_kind real(r_kind),parameter:: four_thirds = 4.0_r_kind / 3.0_r_kind ! Declare local variables logical good,outside,good0,lexist1,lexist2 logical lexist3 character(10) date character(80) hdrstr(2),datstr(2) character(8) subset,subset_check(3) character(30) outmessage character(255) filename integer(i_kind) lnbufr,i,j,k,maxobs,icntpnt,iiout,n,istop integer(i_kind) nmrecs,ibadazm,ibadtilt,ibadrange,ibadwnd,ibaddist,ibadheight,ibadvad,kthin integer(i_kind) iyr,imo,idy,ihr,imn,isc,ithin integer(i_kind) ibadstaheight,ibaderror,notgood,idate,iheightbelowsta,ibadfit integer(i_kind) notgood0 integer(i_kind) novadmatch,ioutofvadrange integer(i_kind) iy,im,idd,ihh,iret,levs,mincy,minobs,kx0,kxadd,kx integer(i_kind) nreal,nchanl,ilat,ilon,ikx integer(i_kind),dimension(5):: idate5 integer(i_kind) ivad,ivadz,nvad,idomsfc integer(i_kind) idbz_bad, idbz_clr, idbz_cld, ibadwnd_cld integer(i_kind) ioutside_rdr, ioutside_obs, ioutside_tim, iazmoutrange integer(i_kind) icld_goodwnd, icld_badwnd, idbz_bad2 integer(i_kind) itiltoutrange, iohgtoutrange, iorngoutrange integer(i_kind) itiltoutrange2 integer(i_kind) nmsg,ntb integer(i_kind) ireadmg, ireadsb integer(i_kind) :: igoodwnd integer(i_kind), dimension(30) :: igoodwnd_dbz real(r_kind) timeb,rmesh,usage,ff10,sfcr,skint,t4dv,t4dvo,toff real(r_kind) eradkm,dlat_earth,dlon_earth real(r_kind) dlat_earth_deg,dlon_earth_deg real(r_kind) dlat,dlon,staheight,tiltangle,clon,slon,clat,slat real(r_kind) timeo,clonh,slonh,clath,slath,cdist,dist real(r_kind) rwnd,azm,height,error,wqm real(r_kind) azm_earth,cosazm_earth,sinazm_earth,cosazm,sinazm real(r_kind):: zsges real(r_kind) dhdr,dsdr,rdbz,range1 real(r_kind) sfcrng0, elvang0 real(r_double),dimension(12):: hdr2 real(r_double),dimension(12):: hdr3 ! ====================================================================! ! variables used for L2RWBUFR data (column-tilt RW) ! ! --- BUFR Key word variables for l2rwbufr_colmtilt by CAPS f88d2arps character(80):: hdrstr2(3), obstr real(r_kind) :: thiserr_rw real(r_double),dimension(5,maxlevs):: radar_obs_columntilt ! to read in column-tilt radar obs in l2rwbufr ! L2 radar data number statistics character(4), dimension(1000) :: radartable_id integer(i_kind), dimension(1000) :: radartable_inside integer(i_kind), dimension(1000) :: radartable_ncolm integer(i_kind), dimension(1000) :: radartable_ncolm_ins integer(i_kind), dimension(1000) :: radartable_ncolm_out integer(i_kind), dimension(1000) :: radartable_ncolm_outtimw integer(i_kind), dimension(1000) :: radartable_ndata_use integer(i_kind) :: icnt_radar_read integer(i_kind) :: icnt_radar_inside integer(i_kind) :: icnt_radar_outside integer(i_kind) :: irdr, irdr_ptr, icnt integer(i_kind) :: ntot_col, ntot_col_o, ntot_col_i, ntot_dat_u, ntot_col_ot integer(i_kind) :: irpt ! filename for level 2 radar obs data on column-tilt (produced by 88d2arps) character(255) :: fname_l2rwbufr_columntilt ! ====================================================================! real(r_kind) thisazm, azm_east integer(i_kind) ncolumntiltl2_in,ncolumntiltl2_kept real(r_kind),dimension(maxdat):: cdata real(r_kind),allocatable,dimension(:,:):: cdata_all real(r_double) rstation_id real(r_double),dimension(12):: hdr character(8) cstaid character(4) this_staid equivalence (this_staid,cstaid) equivalence (cstaid,rstation_id) real(r_double),dimension(7,maxlevs):: radar_obs real(r_double),dimension(4,maxlevs):: vad_obs real(r_double),dimension(2,maxlevs):: fcst_obs character(8) vadid(maxvad) real(r_kind) vadlat(maxvad),vadlon(maxvad),vadqm(maxvad,maxvadbins) real(r_kind) vadu(maxvad,maxvadbins),vadv(maxvad,maxvadbins) real(r_kind) vadcount(maxvad,maxvadbins) real(r_kind),dimension(maxvad,maxvadbins)::vadfit2,vadcount2,vadwgt2 real(r_kind),dimension(maxvad,maxvadbins)::vadfit2_5,vadcount2_5,vadwgt2_5 real(r_kind),dimension(maxvad,maxvadbins)::vadfit3,vadcount3,vadwgt3 real(r_kind) zob,vadqmmin,vadqmmax integer(i_kind) level2(maxvad),level2_5(maxvad),level3(maxvad),level3_tossed_by_2_5(maxvad) integer(i_kind) loop,numcut integer(i_kind) numhits(0:maxvad) real(r_kind) timemax,timemin,errmax,errmin real(r_kind) dlatmax,dlonmax,dlatmin,dlonmin real(r_kind) xscale,xscalei integer(i_kind) max_rrr,nboxmax integer(i_kind) irrr,iaaa,iaaamax,iaaamin integer(i_byte),allocatable::nobs_box(:,:,:,:) real(r_kind) dlonvad,dlatvad,vadlon_earth,vadlat_earth real(r_kind) this_stalat,this_stalon,this_stahgt,thistime,thislat,thislon real(r_kind) azm0,elev0,range0,rotang real(r_kind) thishgt,thisvr,corrected_azimuth,thiserr,corrected_tilt integer(i_kind) nsuper2_in,nsuper2_kept integer(i_kind) nsuper2_5_in,nsuper2_5_kept integer(i_kind) nsuper3_in,nsuper3_kept real(r_kind) errzmax real(r_kind) thisfit,thisvadspd,thisfit2,uob,vob,thiswgt ! real(r_kind) dist2min,dist2max ! real(r_kind) dist2_5min,dist2_5max real(r_kind) vad_leash ! following variables are use for tdr rw data real(r_double),dimension(4,maxlevs):: tdr_obs integer(i_kind) :: ii,jjj,nmissing,nirrr,noutside,ntimeout,nsubzero,iimax integer(i_kind) ntdrvr_in,ntdrvr_kept,ntdrvr_thin1,ntdrvr_thin2 integer(i_kind) ntdrvr_thin2_foreswp,ntdrvr_thin2_aftswp integer(i_kind) maxout,maxdata integer(i_kind) kk,klon1,klat1,klonp1,klatp1 integer(i_kind),allocatable,dimension(:):: isort real(r_single) elevmax,elevmin real(r_single) thisrange,thisazimuth,thistilt real(r_single), dimension(maxlevs) :: dopbin, z, elev, elat8, elon8, glob_azimuth8 real(r_kind) rlon0,this_stalatr,thistiltr real(r_kind) clat0,slat0 real(r_single) a43,aactual,selev0,celev0,erad real(r_kind) sin2,termg,termr,termrg,zobs real(r_kind) xmesh,pmesh real(r_kind),dimension(nsig):: zges,hges real(r_kind) dx,dy,dx1,dy1,w00,w10,w01,w11 logical luse integer(i_kind) ntmp,iout integer(i_kind):: zflag integer(i_kind) nlevz ! vertical level for thinning real(r_kind) crit1,timedif real(r_kind),allocatable,dimension(:):: zl_thin real(r_kind),parameter:: r16000 = 16000.0_r_kind real(r_kind) diffuu,diffvv ! following variables are for fore/aft separation real(r_kind) tdrele1,tdrele2,tdrele3 integer(i_kind) nswp,firstbeam,nforeswp,naftswp,nfore,naft,nswptype,irec logical foreswp,aftswp data lnbufr/10/ data hdrstr(1) / 'CLAT CLON SELV ANEL YEAR MNTH DAYS HOUR MINU MGPT' / data hdrstr(2) / 'PTID YEAR MNTH DAYS HOUR MINU SECO CLAT CLON HSMSL ANAZ ANEL' / data datstr(1) / 'STDM SUPLAT SUPLON HEIT RWND RWAZ RSTD' / data datstr(2) / 'DIST HREF DMVR DVSW' / ! used for BUFR Radial Wind on column-tilt data hdrstr2(1) /'SSTN CLAT CLON HSMSL HSALG ANAZ XOB YOB'/ data hdrstr2(2) /'YEAR MNTH DAYS HOUR MINU SECO'/ data hdrstr2(3) /'SCID HNQV VOCP VOID'/ data obstr /'ANEL HEIT DISTL2RW DMVR HREF'/ data ithin / -9 / data rmesh / -99.999_r_kind / logical print_verbose !*********************************************************************************** print_verbose=.false. if(verbose)print_verbose=.true. ! Check to see if radar wind files exist. If none exist, exit this routine. inquire(file='radar_supobs_from_level2',exist=lexist1) inquire(file=trim(infile),exist=lexist2) lexist3=.false. if ( l_use_rw_columntilt) then ! add level 2 bufr format RW data on column-tilt fname_l2rwbufr_columntilt='l2rwbufr_cltl' inquire(file=trim(fname_l2rwbufr_columntilt),exist=lexist3) if ( .not. lexist3 ) then write(6,'(1x,A70,A20,A10,I6)') & ' level-2 column-tilt radar wind file is not found and to skip : ', & trim(fname_l2rwbufr_columntilt), ' on pe:',mype else write(6,'(1x,A70,A20,A10,I6)') & ' level-2 column-tilt radar wind file is found and to read : ', & trim(fname_l2rwbufr_columntilt), ' on pe:',mype end if end if if (.not.lexist1 .and. .not.lexist2 .and. .not.lexist3) return eradkm=rearth*0.001_r_kind maxobs=2e8 nreal=maxdat nchanl=0 ilon=2 ilat=3 iaaamax=-huge(iaaamax) iaaamin=huge(iaaamin) dlatmax=-huge(dlatmax) dlonmax=-huge(dlonmax) dlatmin=huge(dlatmin) dlonmin=huge(dlonmin) if(ianldate > 2016092000)then hdrstr(2)='PTID YEAR MNTH DAYS HOUR MINU SECO CLAT CLON FLVLST ANAZ ANEL' end if allocate(cdata_all(maxdat,maxobs),isort(maxobs)) isort = 0 cdata_all=zero if (trim(infile) /= 'tldplrbufr' .and. trim(infile) /= 'tldplrso') then ! Initialize variables ! vad_leash=.1_r_kind vad_leash=.3_r_kind ! xscale=5000._r_kind ! xscale=10000._r_kind xscale=20000._r_kind if(print_verbose)then write(6,*)'READ_RADAR: set vad_leash,xscale=',vad_leash,xscale write(6,*)'READ_RADAR: set maxvadbins,maxbadbins*dzvad=',maxvadbins,& maxvadbins*dzvad end if xscalei=one/xscale max_rrr=nint(100000.0_r_kind*xscalei) nboxmax=1 kx0=22500 nmrecs=0 irec=0 errzmax=zero nvad=0 vadlon=zero vadlat=zero vadqm=-99999_r_kind vadu=zero vadv=zero vadcount=zero vadqmmax=-huge(vadqmmax) vadqmmin=huge(vadqmmin) ! First read in all vad winds so can use vad wind quality marks to decide ! which radar data to keep ! Open, then read bufr data open(lnbufr,file=vadfile,form='unformatted') call openbf(lnbufr,'IN',lnbufr) call datelen(10) loop0: do call readsb(lnbufr,iret) if(iret/=0) then call readmg(lnbufr,subset,idate,iret) if(iret/=0) exit loop0 cycle loop0 end if call ufbint(lnbufr,hdr,7,1,levs,'SID XOB YOB DHR TYP SAID TSB') kx=nint(hdr(5)) if(kx /= 224)cycle loop0 ! for now just hardwire vad wind type if(kx==224 .and. .not.newvad) then if(hdr(7)==2) then newvad=.true. exit loop0 end if end if ! End of bufr read loop end do loop0 call closbf(lnbufr) close(lnbufr) ! enddo msg_report open(lnbufr,file=vadfile,form='unformatted') call openbf(lnbufr,'IN',lnbufr) call datelen(10) call readmg(lnbufr,subset,idate,iret) if(iret==0) then ! Time offset call time_4dvar(idate,toff) write(date,'( i10)') idate read (date,'(i4,3i2)') iy,im,idd,ihh if(print_verbose) & write(6,*)'READ_RADAR: first read vad winds--use vad quality marks to qc 2.5/3 radar winds' ! Big loop over vadwnd bufr file loop1: do call readsb(lnbufr,iret) if(iret/=0) then call readmg(lnbufr,subset,idate,iret) if(iret/=0) exit loop1 cycle loop1 end if nmrecs = nmrecs+1 ! Read header. Extract station infomration call ufbint(lnbufr,hdr,7,1,levs,'SID XOB YOB DHR TYP SAID TSB') kx=nint(hdr(5)) if(kx /= 224) cycle loop1 ! for now just hardwire vad wind type ! write(6,*)'new vad::',newvad, hdr(7) if(.not.newvad .and. hdr(7)==2) cycle loop1 if(newvad .and. hdr(7)/=2) cycle loop1 ! and don't worry about subtypes ! Is vadwnd in convinfo file ikx=0 do i=1,nconvtype if(kx == ictype(i)) then ikx=i exit end if end do if(ikx == 0) cycle loop1 ! Time check t4dv=toff+hdr(4) if (l4dvar.or.l4densvar) then if (t4dvwinlen) cycle loop1 ! outside time window else timeb=hdr(4) if(abs(timeb) > ctwind(ikx) .or. abs(timeb) > half) cycle loop1 ! outside time window endif ! Create table of vad lat-lons and quality marks in 500m increments ! for cross-referencing bird qc against radar winds rstation_id=hdr(1) !station id dlon_earth=hdr(2) !station lat (degrees) dlat_earth=hdr(3) !station lon (degrees) if (dlon_earth>=r360) dlon_earth=dlon_earth-r360 if (dlon_earth0) then do i=1,nvad if(modulo(rad2deg*abs(dlon_earth-vadlon(i)),r360)maxvad) then write(6,*)'READ_RADAR: ***ERROR*** MORE THAN ',maxvad,' RADARS: PROGRAM STOPS' call stop2(84) end if ivad=nvad vadlon(ivad)=dlon_earth vadlat(ivad)=dlat_earth vadid(ivad)=cstaid end if ! Update vadqm table call ufbint(lnbufr,vad_obs,4,maxlevs,levs,'ZOB WQM UOB VOB ') call ufbint(lnbufr,fcst_obs,2,maxlevs,levs,'UFC VFC ') if(levs>maxlevs) then write(6,*)'READ_RADAR: ***ERROR*** need to increase read_radar bufr size since ',& ' number of levs=',levs,' > maxlevs=',maxlevs call stop2(84) endif do k=1,levs wqm=vad_obs(2,k) zob=vad_obs(1,k) uob=vad_obs(3,k) vob=vad_obs(4,k) if(newvad) then diffuu=uob-fcst_obs(1,k) diffvv=vob-fcst_obs(2,k) if(sqrt(diffuu**2+diffvv**2)>10.0) cycle if(abs(diffvv)>8.0) cycle if(abs(diffvv)>5.0.and.zob<5000.0) cycle if(zob>7000.0) cycle end if ivadz=nint(zob/dzvad) if(ivadz<1.or.ivadz>maxvadbins) cycle errzmax=max(abs(zob-ivadz*dzvad),errzmax) vadqm(ivad,ivadz)=max(vadqm(ivad,ivadz),wqm) vadqmmax=max(vadqmmax,wqm) vadqmmin=min(vadqmmin,wqm) vadu(ivad,ivadz)=vadu(ivad,ivadz)+uob vadv(ivad,ivadz)=vadv(ivad,ivadz)+vob vadcount(ivad,ivadz)=vadcount(ivad,ivadz)+one end do ! End of bufr read loop end do loop1 ! Normal exit end if call closbf(lnbufr) close(lnbufr) ! Print vadwnd table if(nvad>0) then do ivad=1,nvad do ivadz=1,maxvadbins vadu(ivad,ivadz)=vadu(ivad,ivadz)/max(one,vadcount(ivad,ivadz)) vadv(ivad,ivadz)=vadv(ivad,ivadz)/max(one,vadcount(ivad,ivadz)) end do if(print_verbose) & write(6,'(" n,lat,lon,qm=",i3,2f8.2,2x,25i3)') & ivad,vadlat(ivad)*rad2deg,vadlon(ivad)*rad2deg,(max(-9,nint(vadqm(ivad,k))),k=1,maxvadbins) end do end if if(print_verbose)write(6,*)' errzmax=',errzmax ! Allocate thinning grids around each radar ! space needed is nvad*max_rrr*max_rrr*8*max_zzz ! ! max_rrr=20 ! maxvadbins=20 ! nvad=150 ! space=150*20*20*8*20 = 64000*150=9600000 peanuts allocate(nobs_box(max_rrr,8*max_rrr,maxvadbins,nvad)) nobs_box=0 ! Set level2_5 to 0. Then loop over routine twice, first looking for ! level 2.5 data, and setting level2_5=count of 2.5 data for any 2.5 data ! available that passes the vad tests. The second pass puts in level 3 ! data where it is available and no level 2.5 data was saved/available ! (level2_5=0) vadfit2=zero vadfit2_5=zero vadfit3=zero vadwgt2=zero vadwgt2_5=zero vadwgt3=zero vadcount2=zero vadcount2_5=zero vadcount3=zero level2=0 level2_5=0 level3=0 level3_tossed_by_2_5=0 subset_check(1)='NC006002' subset_check(2)='NC006001' ! First process any level 2 superobs. ! Initialize variables. ikx=0 do i=1,nconvtype if(trim(ioctype(i)) == trim(obstype))ikx = i end do timemax=-huge(timemax) timemin=huge(timemin) errmax=-huge(errmax) errmin=huge(errmin) loop=0 numhits=0 ibadazm=0 ibadwnd=0 ibaddist=0 ibadheight=0 ibadstaheight=0 iheightbelowsta=0 ibaderror=0 ibadvad=0 ibadfit=0 ioutofvadrange=0 kthin=0 novadmatch=0 notgood=0 notgood0=0 nsuper2_in=0 nsuper2_kept=0 ! LEVEL_TWO_READ: if(loop==0 .and. sis=='l2rw') then if(loop==0) outmessage='level 2 superobs:' ! Open sequential file containing superobs open(lnbufr,file='radar_supobs_from_level2',form='unformatted') rewind lnbufr ! dist2max=-huge(dist2max) ! dist2min=huge(dist2min) ! Loop to read superobs data file do read(lnbufr,iostat=iret)this_staid,this_stalat,this_stalon,this_stahgt, & thistime,thislat,thislon,thishgt,thisvr,corrected_azimuth,thiserr,corrected_tilt if(iret/=0) exit nsuper2_in=nsuper2_in+1 dlat_earth=this_stalat !station lat (degrees) dlon_earth=this_stalon !station lon (degrees) if (dlon_earth>=r360) dlon_earth=dlon_earth-r360 if (dlon_earthwinlen) cycle else timeo=thistime if(abs(timeo)>half ) cycle endif ! Get observation (lon,lat). Compute distance from radar. dlat_earth=thislat dlon_earth=thislon if(dlon_earth>=r360) dlon_earth=dlon_earth-r360 if(dlon_earthmax_rrr) cycle ! Extract radial wind data height= thishgt rwnd = thisvr azm_earth = corrected_azimuth if(regional .and. .not. fv3_regional) then cosazm_earth=cos(azm_earth*deg2rad) sinazm_earth=sin(azm_earth*deg2rad) call rotate_wind_ll2xy(cosazm_earth,sinazm_earth,cosazm,sinazm,dlon_earth,dlon,dlat) azm=atan2(sinazm,cosazm)*rad2deg else azm=azm_earth end if iaaa=azm/(r360/(r8*irrr)) iaaa=mod(iaaa,8*irrr) if(iaaa<0) iaaa=iaaa+8*irrr iaaa=iaaa+1 iaaamax=max(iaaamax,iaaa) iaaamin=min(iaaamin,iaaa) error = erradar_inflate*thiserr errmax=max(error,errmax) if(thiserr>zero) errmin=min(error,errmin) ! Perform limited qc based on azimuth angle, radial wind ! speed, distance from radar site, elevation of radar, ! height of observation, observation error, and goodness of fit to vad wind good0=.true. if(abs(azm)>r400) then ibadazm=ibadazm+1; good0=.false. end if if(abs(rwnd)>r200) then ibadwnd=ibadwnd+1; good0=.false. end if if(dist>r400) then ibaddist=ibaddist+1; good0=.false. end if if(staheight<-r1000.or.staheight>r50000) then ibadstaheight=ibadstaheight+1; good0=.false. end if if(height<-r1000.or.height>r50000) then ibadheight=ibadheight+1; good0=.false. end if if(heightr6 .or. thiserr<=zero) then ibaderror=ibaderror+1; good0=.false. end if good=.true. if(.not.good0) then notgood0=notgood0+1 cycle else ! Check fit to vad wind and vad wind quality mark ivadz=nint(thishgt/dzvad) if(ivadz>maxvadbins.or.ivadz<1) then ioutofvadrange=ioutofvadrange+1 cycle end if thiswgt=one/max(r4_r_kind,thiserr**2) thisfit2=(vadu(ivad,ivadz)*cos(azm_earth*deg2rad)+vadv(ivad,ivadz)*sin(azm_earth*deg2rad)-thisvr)**2 thisfit=sqrt(thisfit2) thisvadspd=sqrt(vadu(ivad,ivadz)**2+vadv(ivad,ivadz)**2) vadfit2(ivad,ivadz)=vadfit2(ivad,ivadz)+thiswgt*thisfit2 vadcount2(ivad,ivadz)=vadcount2(ivad,ivadz)+one vadwgt2(ivad,ivadz)=vadwgt2(ivad,ivadz)+thiswgt if(thisfit/max(one,thisvadspd)>vad_leash) then ibadfit=ibadfit+1; good=.false. end if if(nobs_box(irrr,iaaa,ivadz,ivad)>nboxmax) then kthin=kthin+1 good=.false. end if if(vadqm(ivad,ivadz) > r3_5 .or. vadqm(ivad,ivadz) < -one) then ibadvad=ibadvad+1 ; good=.false. end if end if ! If data is good, load into output array if(good) then nsuper2_kept=nsuper2_kept+1 level2(ivad)=level2(ivad)+1 nobs_box(irrr,iaaa,ivadz,ivad)=nobs_box(irrr,iaaa,ivadz,ivad)+1 ndata =min(ndata+1,maxobs) nodata =min(nodata+1,maxobs) !number of obs not used (no meaning here) usage = zero if(icuse(ikx) < 0)usage=r100 if(ncnumgrp(ikx) > 0 )then ! cross validation on if(mod(ndata,ncnumgrp(ikx))== ncgroup(ikx)-1)usage=ncmiter(ikx) end if call deter_sfc2(dlat_earth,dlon_earth,t4dv,idomsfc,skint,ff10,sfcr) LEVEL_TWO_READ: if(loop==0 .and. sis=='l2rw') then cdata(1) = error ! wind obs error (m/s) cdata(2) = dlon ! grid relative longitude cdata(3) = dlat ! grid relative latitude cdata(4) = height ! obs absolute height (m) cdata(5) = rwnd ! wind obs (m/s) cdata(6) = azm*deg2rad ! azimuth angle (radians) cdata(7) = t4dv ! obs time (hour) cdata(8) = ikx ! type cdata(9) = tiltangle ! tilt angle (radians) cdata(10)= staheight ! station elevation (m) cdata(11)= rstation_id ! station id cdata(12)= usage ! usage parameter cdata(13)= idomsfc ! dominate surface type cdata(14)= skint ! skin temperature cdata(15)= ff10 ! 10 meter wind factor cdata(16)= sfcr ! surface roughness cdata(17)=dlon_earth_deg ! earth relative longitude (degrees) cdata(18)=dlat_earth_deg ! earth relative latitude (degrees) cdata(19)=dist ! range from radar in km (used to estimate beam spread) cdata(20)=zsges ! model elevation at radar site cdata(21)=thiserr cdata(22)=two ! if(vadid(ivad)=='0303LWX') then ! dist2max=max(dist2max,dist) ! dist2min=min(dist2min,dist) ! end if do i=1,maxdat cdata_all(i,ndata)=cdata(i) end do END IF LEVEL_TWO_READ else notgood = notgood + 1 end if end do close(lnbufr) ! A simple unformatted fortran file should not be mixed with a bufr I/O LEVEL_TWO_READ_2: if(loop==0 .and. sis=='l2rw') then write(6,*)'READ_RADAR: ',trim(outmessage),' reached eof on 2/2.5/3 superob radar file' write(6,*)'READ_RADAR: nsuper2_in,nsuper2_kept=',nsuper2_in,nsuper2_kept write(6,*)'READ_RADAR: # no vad match =',novadmatch write(6,*)'READ_RADAR: # out of vadrange=',ioutofvadrange write(6,*)'READ_RADAR: # bad azimuths=',ibadazm write(6,*)'READ_RADAR: # bad winds =',ibadwnd write(6,*)'READ_RADAR: # bad dists =',ibaddist write(6,*)'READ_RADAR: # bad stahgts =',ibadstaheight write(6,*)'READ_RADAR: # bad obshgts =',ibadheight write(6,*)'READ_RADAR: # bad errors =',ibaderror write(6,*)'READ_RADAR: # bad vadwnd =',ibadvad write(6,*)'READ_RADAR: # bad fit =',ibadfit write(6,*)'READ_RADAR: # num thinned =',kthin write(6,*)'READ_RADAR: # notgood0 =',notgood0 write(6,*)'READ_RADAR: # notgood =',notgood write(6,*)'READ_RADAR: # hgt belowsta=',iheightbelowsta write(6,*)'READ_RADAR: timemin,max =',timemin,timemax write(6,*)'READ_RADAR: errmin,max =',errmin,errmax write(6,*)'READ_RADAR: dlatmin,max,dlonmin,max=',dlatmin,dlatmax,dlonmin,dlonmax write(6,*)'READ_RADAR: iaaamin,max,8*max_rrr =',iaaamin,iaaamax,8*max_rrr END IF LEVEL_TWO_READ_2 LEVEL_THREE_READ: if(sis=='l3rw' .or. sis=='rw') then ! Next process level 2.5 and 3 superobs ! Bigger loop over first level 2.5 data, and then level3 data timemax=-huge(timemax) timemin=huge(timemin) errmax=-huge(errmax) errmin=huge(errmin) nsuper2_5_in=0 nsuper3_in=0 nsuper2_5_kept=0 nsuper3_kept=0 do loop=1,2 numhits=0 ibadazm=0 ibadwnd=0 ibaddist=0 ibadheight=0 ibadstaheight=0 iheightbelowsta=0 ibaderror=0 ibadvad=0 ibadfit=0 ioutofvadrange=0 kthin=0 novadmatch=0 notgood=0 notgood0=0 ! dist2_5max=-huge(dist2_5max) ! dist2_5min=huge(dist2_5min) if(loop==1) outmessage='level 2.5 superobs:' if(loop==2) outmessage='level 3 superobs:' idate5(1) = iy ! year idate5(2) = im ! month idate5(3) = idd ! day idate5(4) = ihh ! hour idate5(5) = 0 ! minute call w3fs21(idate5,mincy) nmrecs=0 ! Open, then read bufr data open(lnbufr,file=trim(infile),form='unformatted') call openbf(lnbufr,'IN',lnbufr) call datelen(10) call readmg(lnbufr,subset,idate,iret) if(iret==0) then ! Big loop over bufr file loop2: do call readsb(lnbufr,iret) if(iret/=0) then call readmg(lnbufr,subset,idate,iret) if(iret/=0) exit loop2 cycle loop2 end if if(subset/=subset_check(loop)) then call readmg(lnbufr,subset,idate,iret) if(iret/=0) exit loop2 cycle loop2 end if nmrecs = nmrecs+1 ! Read header. Extract station infomration call ufbint(lnbufr,hdr,10,1,levs,hdrstr(1)) ! rstation_id=hdr(1) !station id write(cstaid,'(2i4)')idint(hdr(1)),idint(hdr(2)) if(cstaid(1:1)==' ')cstaid(1:1)='S' dlat_earth=hdr(1) !station lat (degrees) dlon_earth=hdr(2) !station lon (degrees) if (dlon_earth>=r360) dlon_earth=dlon_earth-r360 if (dlon_earth230.0_r_kind .and. & dlat_earth <54.0_r_kind)then cycle loop2 end if end if end if dlat_earth = dlat_earth * deg2rad dlon_earth = dlon_earth * deg2rad if(regional)then call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside) if (outside) cycle loop2 dlatmax=max(dlat,dlatmax) dlonmax=max(dlon,dlonmax) dlatmin=min(dlat,dlatmin) dlonmin=min(dlon,dlonmin) else dlat = dlat_earth dlon = dlon_earth call grdcrd1(dlat,rlats,nlat,1) call grdcrd1(dlon,rlons,nlon,1) endif clon=cos(dlon_earth) slon=sin(dlon_earth) clat=cos(dlat_earth) slat=sin(dlat_earth) staheight=hdr(3) !station elevation tiltangle=hdr(4)*deg2rad ! Find vad wind match ivad=0 do k=1,nvad cdist=sin(vadlat(k))*slat+cos(vadlat(k))*clat* & (sin(vadlon(k))*slon+cos(vadlon(k))*clon) cdist=max(-one,min(cdist,one)) dist=rad2deg*acos(cdist) if(dist < 0.2_r_kind) then ivad=k exit end if end do numhits(ivad)=numhits(ivad)+1 if(ivad==0) then novadmatch=novadmatch+1 cycle loop2 end if vadlon_earth=vadlon(ivad) vadlat_earth=vadlat(ivad) if(regional)then call tll2xy(vadlon_earth,vadlat_earth,dlonvad,dlatvad,outside) if (outside) cycle loop2 dlatmax=max(dlatvad,dlatmax) dlonmax=max(dlonvad,dlonmax) dlatmin=min(dlatvad,dlatmin) dlonmin=min(dlonvad,dlonmin) else dlatvad = vadlat_earth dlonvad = vadlon_earth call grdcrd1(dlatvad,rlats,nlat,1) call grdcrd1(dlonvad,rlons,nlon,1) endif ! Get model terrain at VAD wind location call deter_zsfc_model(dlatvad,dlonvad,zsges) iyr = hdr(5) imo = hdr(6) idy = hdr(7) ihr = hdr(8) imn = hdr(9) idate5(1) = iyr idate5(2) = imo idate5(3) = idy idate5(4) = ihr idate5(5) = imn ikx=0 do i=1,nconvtype if(trim(ioctype(i)) == trim(obstype))ikx = i end do if(ikx==0) cycle loop2 call w3fs21(idate5,minobs) t4dv=real(minobs-iwinbgn,r_kind)*r60inv if (l4dvar.or.l4densvar) then if (t4dvwinlen) cycle loop2 else timeb = real(minobs-mincy,r_kind)*r60inv ! if (abs(timeb)>twind .or. abs(timeb) > ctwind(ikx)) then if (abs(timeb)>half .or. abs(timeb) > ctwind(ikx)) then ! write(6,*)'READ_RADAR: time outside window ',timeb,' skip this obs' cycle loop2 endif endif ! Go through the data levels call ufbint(lnbufr,radar_obs,7,maxlevs,levs,datstr(1)) if(levs>maxlevs) then write(6,*)'READ_RADAR: ***ERROR*** increase read_radar bufr size since ',& 'number of levs=',levs,' > maxlevs=',maxlevs call stop2(84) endif numcut=0 do k=1,levs if(loop==1) nsuper2_5_in=nsuper2_5_in+1 if(loop==2) nsuper3_in=nsuper3_in+1 nread=nread+1 t4dvo=real(minobs+radar_obs(1,k)-iwinbgn,r_kind)*r60inv timemax=max(timemax,t4dvo) timemin=min(timemin,t4dvo) if(loop==2 .and. ivad> 0 .and. level2_5(ivad)/=0) then level3_tossed_by_2_5(ivad)=level3_tossed_by_2_5(ivad)+1 numcut=numcut+1 cycle end if ! Exclude data if it does not fall within time window if (l4dvar.or.l4densvar) then if (t4dvowinlen) cycle timeo=t4dv else timeo=(real(minobs-mincy,r_kind)+real(radar_obs(1,k),r_kind))*r60inv if(abs(timeo)>twind .or. abs(timeo) > ctwind(ikx)) then ! write(6,*)'READ_RADAR: time outside window ',timeo,& ! ' skip obs ',nread,' at lev=',k cycle end if end if ! Get observation (lon,lat). Compute distance from radar. if(radar_obs(3,k)>=r360) radar_obs(3,k)=radar_obs(3,k)-r360 if(radar_obs(3,k)max_rrr) cycle ! Set observation "type" to be function of distance from radar kxadd=nint(dist*one_tenth) kx=kx0+kxadd ! Extract radial wind data height= radar_obs(4,k) rwnd = radar_obs(5,k) azm_earth = r90-radar_obs(6,k) if(regional .and. .not. fv3_regional) then cosazm_earth=cos(azm_earth*deg2rad) sinazm_earth=sin(azm_earth*deg2rad) call rotate_wind_ll2xy(cosazm_earth,sinazm_earth,cosazm,sinazm,dlon_earth,dlon,dlat) azm=atan2(sinazm,cosazm)*rad2deg else azm=azm_earth end if iaaa=azm/(r360/(r8*irrr)) iaaa=mod(iaaa,8*irrr) if(iaaa<0) iaaa=iaaa+8*irrr iaaa=iaaa+1 iaaamax=max(iaaamax,iaaa) iaaamin=min(iaaamin,iaaa) error = erradar_inflate*radar_obs(7,k) ! Increase error for lev2.5 and lev3 if (wrf_nmm_regional.or.nems_nmmb_regional.or.cmaq_regional.or.wrf_mass_regional & .or. fv3_regional ) then if(dlon_earth*rad2deg>230.0_r_kind .and. & dlat_earth*rad2deg <54.0_r_kind)then error = error+r10 end if end if errmax=max(error,errmax) if(radar_obs(7,k)>zero) errmin=min(error,errmin) ! Perform limited qc based on azimuth angle, radial wind ! speed, distance from radar site, elevation of radar, ! height of observation, observation error. good0=.true. if(abs(azm)>r400) then ibadazm=ibadazm+1; good0=.false. end if if(abs(rwnd)>r200) then ibadwnd=ibadwnd+1; good0=.false. end if if(dist>r400) then ibaddist=ibaddist+1; good0=.false. end if if(staheight<-r1000 .or. staheight>r50000) then ibadstaheight=ibadstaheight+1; good0=.false. end if if(height<-r1000 .or. height>r50000) then ibadheight=ibadheight+1; good0=.false. end if if(heightr6 .or. radar_obs(7,k)<=zero) then ibaderror=ibaderror+1; good0=.false. end if good=.true. if(.not.good0) then notgood0=notgood0+1 cycle else ! Check against vad wind quality mark ivadz=nint(height/dzvad) if(ivadz>maxvadbins.or.ivadz<1) then ioutofvadrange=ioutofvadrange+1 cycle end if thiserr = radar_obs(7,k) thiswgt=one/max(r4_r_kind,thiserr**2) thisfit2=(vadu(ivad,ivadz)*cos(azm_earth*deg2rad)+vadv(ivad,ivadz)*sin(azm_earth*deg2rad)-rwnd)**2 thisfit=sqrt(thisfit2) thisvadspd=sqrt(vadu(ivad,ivadz)**2+vadv(ivad,ivadz)**2) if(loop==1) then vadfit2_5(ivad,ivadz)=vadfit2_5(ivad,ivadz)+thiswgt*thisfit2 vadcount2_5(ivad,ivadz)=vadcount2_5(ivad,ivadz)+one vadwgt2_5(ivad,ivadz)=vadwgt2_5(ivad,ivadz)+thiswgt else vadfit3(ivad,ivadz)=vadfit3(ivad,ivadz)+thiswgt*thisfit2 vadcount3(ivad,ivadz)=vadcount3(ivad,ivadz)+one vadwgt3(ivad,ivadz)=vadwgt3(ivad,ivadz)+thiswgt end if if(thisfit/max(one,thisvadspd)>vad_leash) then ibadfit=ibadfit+1; good=.false. end if if(nobs_box(irrr,iaaa,ivadz,ivad)>nboxmax) then kthin=kthin+1 good=.false. end if if(vadqm(ivad,ivadz)>r3_5 .or. vadqm(ivad,ivadz)<-one) then ibadvad=ibadvad+1 ; good=.false. end if end if ! If data is good, load into output array if(good) then if(loop==1.and.ivad>0) then nsuper2_5_kept=nsuper2_5_kept+1 level2_5(ivad)=level2_5(ivad)+1 end if if(loop==2.and.ivad>0) then nsuper3_kept=nsuper3_kept+1 level3(ivad)=level3(ivad)+1 end if nobs_box(irrr,iaaa,ivadz,ivad)=nobs_box(irrr,iaaa,ivadz,ivad)+1 ndata = min(ndata+1,maxobs) nodata = min(nodata+1,maxobs) !number of obs not used (no meaning here) usage = zero if(icuse(ikx) < 0)usage=r100 if(ncnumgrp(ikx) > 0 )then ! cross validation on if(mod(ndata,ncnumgrp(ikx))== ncgroup(ikx)-1)usage=ncmiter(ikx) end if call deter_sfc2(dlat_earth,dlon_earth,t4dv,idomsfc,skint,ff10,sfcr) cdata(1) = error ! wind obs error (m/s) cdata(2) = dlon ! grid relative longitude cdata(3) = dlat ! grid relative latitude cdata(4) = height ! obs absolute height (m) cdata(5) = rwnd ! wind obs (m/s) cdata(6) = azm*deg2rad ! azimuth angle (radians) cdata(7) = t4dvo ! obs time (hour) cdata(8) = ikx ! type cdata(9) = tiltangle ! tilt angle (radians) cdata(10)= staheight ! station elevation (m) cdata(11)= rstation_id ! station id cdata(12)= usage ! usage parameter cdata(13)= idomsfc ! dominate surface type cdata(14)= skint ! skin temperature cdata(15)= ff10 ! 10 meter wind factor cdata(16)= sfcr ! surface roughness cdata(17)=dlon_earth_deg ! earth relative longitude (degrees) cdata(18)=dlat_earth_deg ! earth relative latitude (degrees) cdata(19)=dist ! range from radar in km (used to estimate beam spread) cdata(20)=zsges ! model elevation at radar site cdata(21)=radar_obs(7,k) ! original error from bufr file if(loop==1) then cdata(22)=2.5_r_kind else cdata(22)=three end if do i=1,maxdat cdata_all(i,ndata)=cdata(i) end do else notgood = notgood + 1 end if ! End of k loop over levs end do ! End of bufr read loop end do loop2 end if ! Normal exit ! Close unit to bufr file call closbf(lnbufr) close(lnbufr) write(6,*)'READ_RADAR: ',trim(outmessage),' reached eof on 2.5/3 superob radar file.' if(loop==1) write(6,*)'READ_RADAR: nsuper2_5_in,nsuper2_5_kept=',nsuper2_5_in,nsuper2_5_kept if(loop==2) write(6,*)'READ_RADAR: nsuper3_in,nsuper3_kept=',nsuper3_in,nsuper3_kept write(6,*)'READ_RADAR: # no vad match =',novadmatch write(6,*)'READ_RADAR: # out of vadrange=',ioutofvadrange write(6,*)'READ_RADAR: # bad azimuths=',ibadazm write(6,*)'READ_RADAR: # bad winds =',ibadwnd write(6,*)'READ_RADAR: # bad dists =',ibaddist write(6,*)'READ_RADAR: # bad stahgts =',ibadstaheight write(6,*)'READ_RADAR: # bad obshgts =',ibadheight write(6,*)'READ_RADAR: # bad errors =',ibaderror write(6,*)'READ_RADAR: # bad vadwnd =',ibadvad write(6,*)'READ_RADAR: # bad fit =',ibadfit write(6,*)'READ_RADAR: # num thinned =',kthin write(6,*)'READ_RADAR: # notgood0 =',notgood0 write(6,*)'READ_RADAR: # notgood =',notgood write(6,*)'READ_RADAR: # hgt belowsta=',iheightbelowsta write(6,*)'READ_RADAR: timemin,max =',timemin,timemax write(6,*)'READ_RADAR: errmin,max =',errmin,errmax write(6,*)'READ_RADAR: dlatmin,max,dlonmin,max=',dlatmin,dlatmax,dlonmin,dlonmax write(6,*)'READ_RADAR: iaaamin,max,8*max_rrr =',iaaamin,iaaamax,8*max_rrr end do ! end bigger loop over first level 2.5, then level 3 radar data END IF LEVEL_THREE_READ !! if ( l_use_rw_columntilt) then !=========================================================================================! ! ------------------------------------------------------- ! ! Next process level 2 radar obs (rw+dbz) on column-tilt ! ! produced by 88d2arps of CAPS ! ! Note: ! ! It is assumed that these data have been ! ! processed with good quality control, ! ! so no VAD match is required for these data. ! ! Horizontally they are on model grid, ! ! vertically on radar scan tilt, ! ! so the thinning is different to level 2/2.5/3 ! ! superobs radar data. (thinning is not ready for ! ! this version. ! ! ------------------------------------------------------- ! ! =========================================================================================! ! Bigger loop over first level 2.5 data, and then level3 data timemax=-huge(timemax) timemin=huge(timemin) errmax=-huge(errmax) errmin=huge(errmin) ncolumntiltl2_in=0 ncolumntiltl2_kept=0 loop=-1 thiserr_rw=oe_rw ! user-defined obs error for radial wind in namelist write(6,*)'READ_RADAR: L2RW radar radial wind obs error (user-definedin namelist)=', & thiserr_rw,'(m/s).' numhits=0 ibadazm=0 ibadwnd=0 ibadwnd_cld=0 igoodwnd=0 igoodwnd_dbz(:)=0 ibaddist=0 ibadheight=0 ibadstaheight=0 iheightbelowsta=0 ibaderror=0 ibadvad=0 ibadfit=0 ioutofvadrange=0 idbz_bad=0 idbz_bad2=0 idbz_clr=0 idbz_cld=0 ioutside_rdr=0 ioutside_tim=0 ioutside_obs=0 iazmoutrange=0 itiltoutrange=0 itiltoutrange2=0 iohgtoutrange=0 iorngoutrange=0 icld_goodwnd=0 icld_badwnd=0 kthin=0 novadmatch=0 notgood=0 notgood0=0 radartable_id(:)='' radartable_inside(:)=0 ! 1: radar is inside ! -1: outside ! 0: initial value radartable_ncolm(:)=0 radartable_ncolm_ins(:)=0 radartable_ncolm_out(:)=0 radartable_ncolm_outtimw(:)=0 radartable_ndata_use(:)=0 icnt_radar_read=0 icnt_radar_inside=0 icnt_radar_outside=0 if (l_azm_east1st) then write(6,*)'READ_RADAR: change azimuth to east as 0 before correct it. :',l_azm_east1st else write(6,*)'READ_RADAR: change azimuth to east as 0 after correct it. :',l_azm_east1st end if if (i_correct_tilt ==1) then write(6,*)'READ_RADAR: use GSI code to correct tilt angle. :',i_correct_tilt else write(6,*)'READ_RADAR: use ARPS code to correct tilt angle. :',i_correct_tilt end if if(loop==-1) outmessage='level 2 radar obs (column-tilt):' ! Open, then read bufr data open(lnbufr,file=trim(fname_l2rwbufr_columntilt),form='unformatted') call openbf(lnbufr,'IN',lnbufr) call datelen(10) nmsg=0 ntb=0 nmrecs=0 msg_report: do while (ireadmg(lnbufr,subset,idate) == 0) nmsg=nmsg+1 if(subset(1:6)/='NC0062') then write(6,*)'mype=',mype,' subset check does not pass for L2 columntilt RW file. subset=',subset cycle msg_report end if if (lvldbg > 1) & write(6,'(1x,A6,I4,A30,A10,1x,A12,I10,1x,A10,I8)')'(mype:',mype,') READ_RADAR: subset:', & subset,' data_date:',idate,' nmsg:',nmsg ! Big loop over bufr file sb_report: do while (ireadsb(lnbufr) == 0) nmrecs = nmrecs+1 ! Read header. Extract station infomration hdr(:)=zero call ufbint(lnbufr,hdr,8,1,levs,hdrstr2(1)) ! Read header. Extract time infomration hdr2(:)=zero call ufbint(lnbufr,hdr2,6,1,levs,hdrstr2(2)) ! Read header. Extract other infomration (not used by GSI currently) hdr3(:)=zero call ufbint(lnbufr,hdr3,4,1,levs,hdrstr2(3)) rstation_id=hdr(1) !station id SSTN if(abs(hdr(2))>r90 .or. abs(hdr(3))>r360) cycle sb_report dlat_earth=hdr(2) !station lon (degrees) CLAT dlon_earth=hdr(3) !station lat (degrees) CLON if (dlon_earth==r360) dlon_earth=dlon_earth-r360 if (dlon_earth 0 ) then write(6,*)'*** WARNING --> READ_RADAR: (columntilt) radar:',this_staid, & ' is found both inside and outside. Wrong! <-- WARNING ***' call stop2(999) end if radartable_inside(irdr_ptr) = -1 radartable_ncolm_out(irdr_ptr) = radartable_ncolm_out(irdr_ptr) + 1 ioutside_obs=ioutside_obs+1 cycle sb_report else if ( radartable_inside(irdr_ptr) < 0 ) then write(6,*)'*** WARNING --> READ_RADAR: (columntilt) radar:',this_staid, & ' is found both inside and outside. Wrong! <-- WARNING ***' call stop2(999) end if radartable_inside(irdr_ptr) = 1 end if dlatmax=max(dlat,dlatmax) dlonmax=max(dlon,dlonmax) dlatmin=min(dlat,dlatmin) dlonmin=min(dlon,dlonmin) else dlat = dlat_earth dlon = dlon_earth call grdcrd1(dlat,rlats,nlat,1) call grdcrd1(dlon,rlons,nlon,1) radartable_inside(irdr_ptr) = 1 ! global: radar is always inside endif clon=cos(dlon_earth) slon=sin(dlon_earth) clat=cos(dlat_earth) slat=sin(dlat_earth) staheight=hdr(4) ! station elevation HSMSL ! ---------------------------------------------------------------------------------! ! Skipping finding vad wind match for column-tilt radar wind data ! ---------------------------------------------------------------------------------! ! Get model terrain at radar station location: zsges call deter_zsfc_model(dlat,dlon,zsges) iyr = hdr2(1) imo = hdr2(2) idy = hdr2(3) ihr = hdr2(4) imn = hdr2(5) isc = hdr2(6) idate5(1) = iyr idate5(2) = imo idate5(3) = idy idate5(4) = ihr idate5(5) = imn ikx=0 do i=1,nconvtype if(trim(ioctype(i)) == trim(obstype))ikx = i end do if(ikx==0) cycle sb_report ! time window check call w3fs21(idate5,minobs) call w3fs21(idate5,mincy) t4dv=real(minobs-iwinbgn,r_kind)*r60inv if (l4dvar.or.l4densvar) then if (t4dvwinlen) cycle sb_report else timeb = real(minobs-mincy,r_kind)*r60inv if (abs(timeb)>half .or. abs(timeb) > ctwind(ikx)) then ioutside_tim=ioutside_tim+1 radartable_ncolm_outtimw(irdr_ptr) = radartable_ncolm_outtimw(irdr_ptr) + 1 cycle sb_report endif endif ! Get observation (lon,lat). Compute distance from radar. ! adjusting XOB (longitude) --> hdr(7) ! YOB (lattitude) --> hdr(8) if(abs(hdr(8))>r90 .or. abs(hdr(7))>r360) cycle sb_report thislon=hdr(7) ! obs longitude (Deg) thislat=hdr(8) ! obs latitude (Deg) if(thislon==r360) thislon=thislon-r360 if(thislonr89_5) cycle sb_report ! Checking the azimuth angle and counting the bad data if (hdr(6) < 0.0_r_kind .or. hdr(6) > 360.0_r_kind) then iazmoutrange=iazmoutrange+1 write(6,'(1x,A100,F10.3,1x,I8)')'(mype:',mype, & ') READ_RADAR: azimuth in BUFR is not in 0~360 (skip) -- azm_hdr(6) cnt:', & hdr(6),iazmoutrange cycle sb_report end if ! azm_east is only used for cdata(13) when diagnose and plotting azm_east = r90 - hdr(6) ! original azimuth (degree, but east as 0) if ( azm_east < zero ) azm_east = azm_east + r360 if ( abs(azm_east - r360) <= 1.0E-3_r_kind ) azm_east = zero if ( azm_east > r360 ) azm_east = azm_east - r360 ! Adjustment of Azimuth. ! Including: ! 1) convert azimuth from north as zero (radar-based) to east as zero (model-based) ! 2) correcting azimuth due to the differences between the earth-lat/lon ! of radar station and radar observation ! 3) rotation of azimuth from earth-lat/lon model grid to WRF Curvilinear grid if ( l_azm_east1st ) then ! 1. First, convert azimuth from pointing to north defined as 0 to pointing to east as 0 (GSI rule) ! (1) hdr(6) -- radar azimuth, north as 0 (degree) ! (2) thisazm -- model wind direction azimuth with east as 0 (degree) thisazm = r90-hdr(6) ! wind direction azimuth (deg) ! (3) change the domain of thisazm back to [0,360]. if ( thisazm < zero ) thisazm = thisazm + r360 if ( abs(thisazm - r360) <= 1.0E-3_r_kind ) thisazm = zero if ( thisazm > r360 ) thisazm = thisazm - r360 ! 2. then azimuth is corrected due to lat/lon differences between radar stn and radar obs. ! see subroutine radar_bufr_read_all (in read_l2bufr_mod.f90) for "corrected_azimuth". ! azm_earth is the corrected azimuth (degree) if ( l_correct_azmu ) then call get_azimuth_corrected(this_stalon,this_stalat, & thislon, thislat, & thisazm, azm_earth) ! Note: the corrected_azimuth is in [-pi,pi] (due to domain of Function ATAN2). ! Convert it back to [0,2*pi], i.e. [0, 360] if ( azm_earth < zero) azm_earth = azm_earth + r360 if ( abs(azm_earth - r360) <= 1.0E-3_r_kind ) azm_earth = zero if ( azm_earth > r360 ) azm_earth = azm_earth - r360 else azm_earth = thisazm ! if do not correct it. end if else ! 1. First, azimuth is corrected due to lat/lon differences between radar stn and radar obs. thisazm = hdr(6) ! ANAZ (degree and zero to north) ! see subroutine radar_bufr_read_all (in read_l2bufr_mod.f90) for "corrected_azimuth". ! azm_earth is the corrected azimuth (degree) if ( l_correct_azmu ) then call get_azimuth_corrected(this_stalon,this_stalat, & thislon, thislat, & thisazm, azm_earth) ! Note: the corrected_azimuth is in [-pi,pi] (due to domain of Function ATAN2). ! Convert it back to [0,2*pi], i.e. [0, 360] if ( azm_earth < zero) azm_earth = azm_earth + r360 if ( abs(azm_earth - r360) <= 1.0E-3_r_kind ) azm_earth = zero if ( azm_earth > r360 ) azm_earth = azm_earth - r360 else azm_earth = thisazm ! if do not correct it. end if ! 2. Then convert azimuth from pointing to north defined as 0 to pointing to east as 0 (GSI rule) azm_earth = r90-azm_earth ! east as zero, (degree) ! (1) change the domain of azm_earth back to [0,360]. if ( azm_earth < zero ) azm_earth = azm_earth + r360 if ( abs(azm_earth - r360) <= 1.0E-3_r_kind ) azm_earth = zero if ( azm_earth > r360 ) azm_earth = azm_earth - r360 end if ! 3. rotation adjustment due to lat/lon earth coordinate rotating to X/Y ! rotated coordinate. So rotaing azimuth_eaarth to azimuth_xy ! Note: ! (1) for this calculation, azimuth must be changed to the angle ! with 0 to the east, not to the north. (r90-azm_earth) ! (2) after rotation, azm (the rotated azimuth) is in unit of ! degree. So when stored in cdata(6), it is converted to ! radians. if(regional .and. .not. fv3_regional) then cosazm_earth=cos(azm_earth*deg2rad) sinazm_earth=sin(azm_earth*deg2rad) call rotate_wind_ll2xy(cosazm_earth,sinazm_earth,cosazm,sinazm,dlon_earth,dlon,dlat) azm=atan2(sinazm,cosazm)*rad2deg ! Note: ATAN2 is used to find the rotated_azimuth, the domain of ! ATAN2 is [-pi,pi], or say [-180, 180]. For consistency, ! convert it back to [0,2*pi],i.e., [0, 360]. if (azm < zero) azm = azm + r360 if ( abs(azm - r360) <= 1.0E-3_r_kind ) azm = zero if ( azm > r360 ) azm = azm - r360 else azm=azm_earth end if ! Go through the data levels call ufbint(lnbufr,radar_obs_columntilt,5,maxlevs,levs,obstr) if(levs>maxlevs) then write(6,*)'READ_RADAR: ***ERROR*** increase read_radar bufr size since ',& 'number of levs=',levs,' > maxlevs=',maxlevs call stop2(84) endif numcut=0 ! l3 data toosed by l2.5 data loop_lvl : do k=1,levs ncolumntiltl2_in=ncolumntiltl2_in+1 nread=nread+1 ntb=ntb+1 t4dvo=real(minobs-iwinbgn,r_kind)*r60inv timemax=max(timemax,t4dvo) timemin=min(timemin,t4dvo) kxadd=nint(dist*one_tenth) kx=kx0+kxadd ! Extract radial wind data thistilt = radar_obs_columntilt(1,k) ! radar scan/tilt angle (degree) height = radar_obs_columntilt(2,k) ! obs beam height (meter) range1 = radar_obs_columntilt(3,k) ! obs beam range (meter) rwnd = radar_obs_columntilt(4,k) ! radial wind (m/s) rdbz = radar_obs_columntilt(5,k) ! reflectivity (dBZ) sfcrng0 = dist*r1000 ! surface range in unit of meters ! check on the values from BUFR data ! Should check dbz and rwnd also (first?) if (thistilt < 0.0_r_kind .OR. thistilt > 90.0_r_kind) then itiltoutrange=itiltoutrange+1 if ( rdbz >= 0.0_r_kind .AND. rdbz < 100.0_r_kind ) then write(9997,'(1x,A80,F10.3,1x,I8,2F20.3)') & 'READ_RADAR: tilt in BUFR is not in 0~90 (skip) -- tilt cnt:', & thistilt,itiltoutrange,rdbz,rwnd end if ! cycle loop_lvl end if if (height < -1.0_r_kind .OR. height > 100000.0_r_kind ) then iohgtoutrange=iohgtoutrange+1 if ( rdbz >= 0.0_r_kind .AND. rdbz < 100.0_r_kind ) then write(9997,'(1x,A80,F20.3,1x,I8,2F20.3)') & 'READ_RADAR: obs hgt in BUFR is not in 0~100km (skip) -- obshgt cnt:', & height,iohgtoutrange,rdbz,rwnd end if ! cycle loop_lvl end if if (range1 < 0.0_r_kind .OR. range1 > 510000.0_r_kind ) then iorngoutrange=iorngoutrange+1 if ( rdbz >= 0.0_r_kind .AND. rdbz < 100.0_r_kind ) then write(9997,'(1x,A80,F20.3,1x,I8,2F20.3)') & 'READ_RADAR: obs range in BUFR is not in 0~510km (skip) -- obsrange cnt:', & range1,iorngoutrange,rdbz,rwnd end if ! cycle loop_lvl end if ! Get corrected tilt(elevation) angle if ( l_correct_tilt ) then if ( i_correct_tilt == 1 ) then ! Considering the curvature effect of beam and earth surface ! Note: two schemes ! 1. algorithm used in radar_bufr_read_all (read_l2bufr_mod.f90) for superobs call get_rdr_obshgttilt(thistilt,range1,staheight, & thishgt,corrected_tilt) write(9991,'(1x,A4,A4,A9,F7.3,A7,F12.1,A9,F7.1,A12,F7.1,A12,F7.1,A11,F7.3,A8,F12.1)') & 'ID:',this_staid,' tilt(o):',thistilt,' range:',range1, & ' sta_hgt:',staheight,' obs_hgt(c):',thishgt,' obs_hgt(o):',height, & ' tilt(crc):',corrected_tilt,' sfcrng:',sfcrng0 else ! 2. scheme used in arps pacakge (arpsenkf: dhdr and dsdr) ! note: ! using beamelv and dhdrange in adas/radarlib3d.f90. Need ! obs height (beam height) and surface range (from radar ! to obs location on surface). Based on 88d2arps data, ! when rdbz/rwnd are missing values, the obs height is ! also missing values (-999.0), then it is not necessary ! to correct the tilt, since this obs will be rejected. call beamelv(height,sfcrng0,elvang0,range0) call dhdrange(elvang0,range0,dhdr) dsdr=sqrt(1._r_kind-dhdr*dhdr) corrected_tilt=atan2(dhdr,dsdr)*rad2deg end if ! Warning Note: tiltangle which is stored in cdata(9), must be in unit of radians. tiltangle=corrected_tilt else tiltangle=thistilt end if if (tiltangle < 0.0_r_kind .OR. tiltangle > 90.0_r_kind) then itiltoutrange2=itiltoutrange2+1 write(9997,'(1x,A100,2F10.3,1x,I8,4F20.3)') & 'READ_RADAR: crc-tilt is not in 0~90deg. -- tilt (org crc) cnt dbz rw height range:', & thistilt, tiltangle, itiltoutrange2, rdbz, rwnd, height, range1 ! cycle lvl_loop end if ! obs error inflation with erradar_inflate (not really used for L2RWBUFR data) error = erradar_inflate*oe_rw ! Increase error for lev2.5 and lev3 ! (for specific domain: but commented off and not used for L2RWBUFR data) errmax=max(error,errmax) if(oe_rw > zero) errmin=min(error,errmin) ! Perform limited qc based on azimuth angle, radial wind ! speed, distance from radar site, elevation of radar, ! height of observation, observation error. good0=.true. if(abs(azm)>r400) then ibadazm=ibadazm+1; good0=.false. end if if(abs(rwnd)>r200) then ibadwnd=ibadwnd+1; good0=.false. if (rdbz >= 5.0_r_kind .and. rdbz < 120._r_kind) then ibadwnd_cld=ibadwnd_cld + 1 write(9998,'(1x,A60,2F20.2,A8,2F10.3,2F10.2,I10)') & 'bad wind with cloud --- rw dbz rdr_id azm scntlt xob yob cnt:', & rwnd,rdbz,this_staid,azm,thistilt,dlon,dlat,ibadwnd_cld end if else igoodwnd=igoodwnd+1 if (rdbz<-20.0_r_kind ) then igoodwnd_dbz(1) = igoodwnd_dbz(1)+1 else if (rdbz>=-20.0_r_kind .AND. rdbz < -0.5_r_kind ) then igoodwnd_dbz(2) = igoodwnd_dbz(2)+1 else if (rdbz>=-0.5_r_kind .AND. rdbz < 0.5_r_kind ) then igoodwnd_dbz(3) = igoodwnd_dbz(3)+1 else if (rdbz>= 0.5_r_kind .AND. rdbz < 5.0_r_kind ) then igoodwnd_dbz(4) = igoodwnd_dbz(4)+1 else if (rdbz>= 5.0_r_kind .AND. rdbz < 100.0_r_kind ) then igoodwnd_dbz(5) = igoodwnd_dbz(5)+1 else igoodwnd_dbz(6) = igoodwnd_dbz(6)+1 end if end if if(dist>r400) then ibaddist=ibaddist+1; good0=.false. end if if(staheight<-r1000 .or. staheight>r50000) then ibadstaheight=ibadstaheight+1; good0=.false. end if if(height<-r1000 .or. height>r50000) then ibadheight=ibadheight+1; good0=.false. end if if(heightr6 .or. oe_rw<=zero) then ibaderror=ibaderror+1; good0=.false. end if if (rdbz < -20.0_r_kind ) then idbz_bad=idbz_bad+1; good0=.false. else if (rdbz >= -20.0_r_kind .AND. rdbz < refl_lowbnd_rw) then idbz_clr=idbz_clr+1; good0=.false. else if (rdbz >= refl_lowbnd_rw .and. rdbz < 100.0_r_kind) then idbz_cld=idbz_cld+1 if (abs(rwnd)>r200) then icld_badwnd=icld_badwnd+1 good0=.false. else icld_goodwnd=icld_goodwnd+1 end if else idbz_bad2=idbz_bad2+1; good0=.false. end if ! good and notgood are used for further QC check with VAD match. ! note: No further QC check with VAD match for L2RWBUFR data ! So good is actually same as good0 good=.true. if(.not.good0) then notgood0=notgood0+1 good=.false. cycle loop_lvl end if ! If data is good, load into output array if(good) then ! no vad and obs box counting for L2RWBUFR ! (they are used for NCEP L2/2.5/3 radar obs, see code for those data) radartable_ndata_use(irdr_ptr) = radartable_ndata_use(irdr_ptr) + 1 ncolumntiltl2_kept=ncolumntiltl2_kept+1 if ( ndata+1 > maxobs ) then do irpt=1,5 write(6,'(1x,A128)')' *** *** WARNING --> READ_RADAR: ndata is more than maxobs(=2,000,000) & --> increase maxobs!!! <-- WARNING *** ***' end do write(6,'(1x,A128)')' *** *** WARNING --> INCREASE maxobs in READ_RADAR beyond 2,000,000, & re-compile GSI, re-run !!! <-- WARNING*** ***' end if ndata = min(ndata+1,maxobs) nodata = min(nodata+1,maxobs) !number of obs not used (no meaning here) usage = zero if(icuse(ikx) < 0)usage=r100 call deter_sfc2(dlat_earth,dlon_earth,t4dv,idomsfc,skint,ff10,sfcr) cdata(1) = error ! wind obs error (m/s) cdata(2) = dlon ! grid relative longitude cdata(3) = dlat ! grid relative latitude cdata(4) = height ! obs absolute height (m) cdata(5) = rwnd ! wind obs (m/s) cdata(6) = azm*deg2rad ! azimuth angle (radians) cdata(7) = t4dvo ! obs time (hour) cdata(8) = ikx ! type cdata(9) = tiltangle*deg2rad ! tilt angle (radians) cdata(10)= staheight ! station elevation (m) cdata(11)= rstation_id ! station id cdata(12)= usage ! usage parameter cdata(13)= idomsfc ! dominate surface type cdata(14)= skint ! skin temperature cdata(15)= ff10 ! 10 meter wind factor cdata(16)= sfcr ! surface roughness cdata(17)=dlon_earth*rad2deg ! earth relative longitude (degrees) cdata(18)=dlat_earth*rad2deg ! earth relative latitude (degrees) cdata(19)=dist ! range from radar in km (used to estimate beam spread) cdata(20)=zsges ! model elevation at radar site cdata(21)=thiserr_rw ! original error (same as user-defined in namelist) cdata(22)=two ! sub-type, rank of radar date (level 2/2.5/3. or higher for airborne) ! put original azimuth, tilt and dbz into cdata for plotting diag file. if (l_plt_diag_rw) then cdata(13)= azm_east ! original azimuth (degree, but east as 0) cdata(14)= thistilt ! orignal scan tilt (degree) cdata(16)= rdbz ! reflectivity at same point with rwnd end if do i=1,maxdat cdata_all(i,ndata)=cdata(i) end do else notgood = notgood + 1 end if ! End of k loop over levs end do loop_lvl ! End of bufr read loop end do sb_report end do msg_report ! Close unit to bufr file call closbf(lnbufr) write(6,*)'READ_RADAR: ',trim(outmessage),' reached eof on level 2 radar file (grid-tilt; column-tilt).' if(loop==-1) write(6,*)'READ_RADAR: ncolumntiltl2_in,ncolumntiltl2_kept=',ncolumntiltl2_in,ncolumntiltl2_kept write(6,*)'READ_RADAR: # no vad match =',novadmatch write(6,*)'READ_RADAR: # out of vadrange=',ioutofvadrange write(6,*)'READ_RADAR: # bad azimuths=',ibadazm write(6,*)'READ_RADAR: # bad winds =',ibadwnd write(6,*)'READ_RADAR: # bad dists =',ibaddist write(6,*)'READ_RADAR: # bad stahgts =',ibadstaheight write(6,*)'READ_RADAR: # bad obshgts =',ibadheight write(6,*)'READ_RADAR: # bad errors =',ibaderror write(6,*)'READ_RADAR: # bad vadwnd =',ibadvad write(6,*)'READ_RADAR: # bad fit =',ibadfit write(6,*)'READ_RADAR: # num thinned =',kthin write(6,*)'READ_RADAR: # hgt belowsta=',iheightbelowsta write(6,*)'READ_RADAR: # bad dbz(<-20.0) =',idbz_bad write(6,*)'READ_RADAR: # bad dbz(>100.0) =',idbz_bad2 write(6,*)'READ_RADAR: # clearsky dbz (-20.0~',refl_lowbnd_rw,') =',idbz_clr write(6,*)'READ_RADAR: # cloud dbz (',refl_lowbnd_rw,'~100.0) =',idbz_cld write(6,*)'READ_RADAR: # bad winds (cloud) =',ibadwnd_cld write(6,*)'READ_RADAR: # cloud (bad winds) =',icld_badwnd write(6,*)'READ_RADAR: # cloud (goodwinds) =',icld_goodwnd write(6,*)'READ_RADAR: # azim out of range =',iazmoutrange write(6,*)'READ_RADAR: # tilt(org) out of range =',itiltoutrange write(6,*)'READ_RADAR: # tilt(crc) out of range =',itiltoutrange2 write(6,*)'READ_RADAR: # obshgt out range =',iohgtoutrange write(6,*)'READ_RADAR: # obsrng out range =',iorngoutrange write(6,*)'READ_RADAR: # notgood0 =',notgood0 write(6,*)'READ_RADAR: # notgood(VAD)=',notgood write(6,*)'READ_RADAR: timemin,max =',timemin,timemax write(6,*)'READ_RADAR: errmin,max =',errmin,errmax write(6,*)'READ_RADAR: dlatmin,max,dlonmin,max=',dlatmin,dlatmax,dlonmin,dlonmax write(6,*)'READ_RADAR: nmsg,nmrecs,ntb=',nmsg,nmrecs,ntb write(6,*)'READ_RADAR: # good winds (abs(rwnd) <= 200m/s) =',igoodwnd write(6,*)'READ_RADAR: # dbz distribution in good winds #' write(6,'(1x,A60,F12.5,1x,F12.5)')'READ_RADAR:radar wind--> obs error(adjust) thiserr_rw(original)= ',error, thiserr_rw write(6,*)'----------------------------------------------------------------' write(6,'(1x,6A15)')'<-20.0dbz','-20.0 ~ -0.5','-0.5 ~ 0.5','0.5 ~ 5.0','5.0 ~ 100.0','>100.0dbz' write(6,'(1x,6I15)')igoodwnd_dbz(1:6) write(6,*)'----------------------------------------------------------------' ! sum up radar table information ntot_col = 0; ntot_col_o = 0; ntot_col_i = 0; ntot_dat_u = 0; ntot_col_ot= 0; icnt_radar_inside=0; icnt_radar_outside=0; open(2018,file='radarlist_rw_used.txt',form='formatted') write(2018,'(1x,7(1x,A15))')' RADAR_ID', ' radar_inside', ' total_column', & ' column_inside', ' column_outside', ' OutTWindow', & ' ndata_use' write(6,*)'----------------------------------------------------------------' write(6, '(1x,7(1x,A15))')' RADAR_ID', ' radar_inside', ' total_column', & ' column_inside', ' column_outside', 'colm_OutTWindow', & ' ndata_use' do icnt=1, icnt_radar_read if ( radartable_inside(icnt) > 0 ) then icnt_radar_inside = icnt_radar_inside + 1 else icnt_radar_outside = icnt_radar_outside + 1 end if ntot_col = ntot_col + radartable_ncolm(icnt) ntot_col_o = ntot_col_o + radartable_ncolm_out(icnt) ntot_col_i = ntot_col_i + radartable_ncolm_ins(icnt) ntot_dat_u = ntot_dat_u + radartable_ndata_use(icnt) ntot_col_ot= ntot_col_ot+ radartable_ncolm_outtimw(icnt) write(2018,'(1x,1x,A15,1x,I15,5(1x,I15))')radartable_id(icnt),radartable_inside(icnt), & radartable_ncolm(icnt),radartable_ncolm_ins(icnt),radartable_ncolm_out(icnt), & radartable_ncolm_outtimw(icnt),radartable_ndata_use(icnt) write(6,'(1x,1x,A15,1x,I15,5(1x,I15))')radartable_id(icnt),radartable_inside(icnt), & radartable_ncolm(icnt),radartable_ncolm_ins(icnt),radartable_ncolm_out(icnt), & radartable_ncolm_outtimw(icnt),radartable_ndata_use(icnt) end do write(2018,*)'----------------------------------------------------------------' write(6,*)'----------------------------------------------------------------' write(2018,'(1x,1x,A15,1x,I6,3x,I6,5(1x,I15))')' Total ',icnt_radar_inside, & icnt_radar_outside,ntot_col,ntot_col_i,ntot_col_o,ntot_col_ot,ntot_dat_u write(6, '(1x,1x,A15,1x,I6,3x,I6,5(1x,I15))')' Total ',icnt_radar_inside, & icnt_radar_outside,ntot_col,ntot_col_i,ntot_col_o,ntot_col_ot,ntot_dat_u write(2018,*)'----------------------------------------------------------------' close(2018) write(6,*)'----------------------------------------------------------------' ! end of read columntilt radar level 2 data !====================================================================================! end if ! Write out vad statistics do ivad=1,nvad if(print_verbose)write(6,'(" fit of 2, 2.5, 3 data to vad station, lat, lon = ",a8,2f14.2)') & vadid(ivad),vadlat(ivad)*rad2deg,vadlon(ivad)*rad2deg do ivadz=1,maxvadbins if(vadcount2(ivad,ivadz) > half .and. vadcount2_5(ivad,ivadz) > half & .and. vadcount(ivad,ivadz) > half)then if(vadcount2(ivad,ivadz)>half) then vadfit2(ivad,ivadz)=sqrt(vadfit2(ivad,ivadz)/vadwgt2(ivad,ivadz)) else vadfit2(ivad,ivadz)=zero end if if(vadcount2_5(ivad,ivadz)>half) then vadfit2_5(ivad,ivadz)=sqrt(vadfit2_5(ivad,ivadz)/vadwgt2_5(ivad,ivadz)) else vadfit2_5(ivad,ivadz)=zero end if if(vadcount3(ivad,ivadz)>half) then vadfit3(ivad,ivadz)=sqrt(vadfit3(ivad,ivadz)/vadwgt3(ivad,ivadz)) else vadfit3(ivad,ivadz)=zero end if if(print_verbose)write(6,'(" h,f2,f2.5,f3=",i7,f10.2,"/",i5,f10.2,"/",i5,f10.2,"/",i5)')nint(ivadz*dzvad),& vadfit2(ivad,ivadz),nint(vadcount2(ivad,ivadz)),& vadfit2_5(ivad,ivadz),nint(vadcount2_5(ivad,ivadz)),& vadfit3(ivad,ivadz),nint(vadcount3(ivad,ivadz)) end if end do end do deallocate(nobs_box) end if erad = rearth thiserr=5.0_r_kind timemax=-huge(timemax) timemin=huge(timemin) errmax=-huge(errmax) errmin=huge(errmin) elevmax=-huge(elevmax) elevmin=huge(elevmin) loop=3 nirrr=0 noutside=0 ntimeout=0 nsubzero=0 ibadazm=0 ibadwnd=0 ibaddist=0 ibadtilt=0 ibadrange=0 ibadheight=0 ibadstaheight=0 notgood=0 notgood0=0 nread=0 ntdrvr_in=0 ntdrvr_kept=0 ntdrvr_thin1=0 ntdrvr_thin2=0 ntdrvr_thin2_foreswp=0 ntdrvr_thin2_aftswp=0 maxout=0 maxdata=0 nmissing=0 subset_check(3)='NC006070' icntpnt=0 nswp=0 nforeswp=0 naftswp=0 nfore=0 naft=0 xscale=100._r_kind xscalei=one/xscale max_rrr=nint(100000.0_r_kind*xscalei) jjj=0 iimax=0 if(loop == 3) outmessage='tail Doppler radar obs:' use_all = .true. do i=1,nconvtype if(trim(ioctype(i)) == trim(obstype) .and. ictype(i) < 999 .and. icuse(i) > 0)then ithin=ithin_conv(i) if(ithin > 0)then rmesh=rmesh_conv(i) pmesh=pmesh_conv(i) use_all = .false. if(pmesh > zero) then ! Here pmesh is height in meters zflag=1 nlevz=r16000/pmesh else zflag=0 nlevz=nsig endif xmesh=rmesh call make3grids(xmesh,nlevz) allocate(zl_thin(nlevz)) if (zflag==1) then do k=1,nlevz zl_thin(k)=(k-1)*pmesh enddo endif write(6,*)'READ_RADAR: obstype,ictype,rmesh,zflag,nlevz,pmesh=',& trim(ioctype(i)),ictype(i),rmesh,zflag,nlevz,pmesh exit end if end if end do if(trim(infile) == 'tldplrso') then ! Loop to read TDR superobs data ikx=0 do i=1,nconvtype if(trim(ioctype(i)) == trim(obstype))ikx = i end do if(ikx == 0) return call w3fs21(iadate,mincy) ! analysis time in minutes open(lnbufr,file=trim(infile),form='formatted',err=300) rewind (lnbufr) do n=1,10 istop=0 read(lnbufr,'(a)',err=200,end=1200)filename print *,'filename=', trim(filename) open(25,file=trim(filename),form='formatted',access='sequential') loop3: do while (istop.eq.0) ii=1 READ(25,'(I4,4I2,8F10.3)',iostat=istop) iyr,imo,idy,ihr,imn,this_stalat, & this_stalon,this_stahgt,azm0,elev0,range0,thisvr,rotang nread=nread+1 idate5(1) = iyr idate5(2) = imo idate5(3) = idy idate5(4) = ihr idate5(5) = imn call w3fs21(idate5,minobs) t4dv=real(minobs-iwinbgn,r_kind)*r60inv if (l4dvar.or.l4densvar) then if (t4dvwinlen) cycle loop3 timeo=t4dv else timeo = real(minobs-mincy,r_kind)*r60inv if (abs(timeo)>twind) cycle loop3 endif timemax=max(timemax,timeo) timemin=min(timemin,timeo) rlon0=deg2rad*this_stalon this_stalatr=this_stalat*deg2rad clat0=cos(this_stalatr) ; slat0=sin(this_stalatr) thistilt=elev0 elevmax=max(elevmax,thistilt) elevmin=min(elevmin,thistilt) thisazimuth=azm0 thisrange=range0*r1000 if(abs(thistilt)>r75)then ibadtilt=ibadtilt+1; cycle loop3 endif staheight=this_stahgt if(staheight<-r1000.or.staheight>r50000) then ibadstaheight=ibadstaheight+1; cycle loop3 end if aactual=erad+this_stahgt thistiltr=thistilt*deg2rad selev0=sin(thistiltr) ; celev0=cos(thistiltr) a43=four_thirds*aactual call getvrlocalinfo(thisrange,thisazimuth,this_stahgt,aactual,a43,selev0,celev0, & rlon0,clat0,slat0,r8,r89_5,nsubzero,ii,z(ii),elev(ii),elat8(ii), & elon8(ii),glob_azimuth8(ii)) dlat_earth=this_stalat !station lat (degrees) dlon_earth=this_stalon !station lon (degrees) if (dlon_earth>=r360) dlon_earth=dlon_earth-r360 if (dlon_earth=r360) dlon_earth=dlon_earth-r360 if(dlon_earthmax_rrr)then nirrr=nirrr+1 cycle endif ! Extract radial wind data height= z(ii) rwnd = thisvr azm_earth = glob_azimuth8(ii) if(regional .and. .not. fv3_regional) then cosazm_earth=cos(azm_earth*deg2rad) sinazm_earth=sin(azm_earth*deg2rad) call rotate_wind_ll2xy(cosazm_earth,sinazm_earth,cosazm,sinazm,dlon_earth,dlon,dlat) azm=atan2(sinazm,cosazm)*rad2deg else azm=azm_earth end if iaaa=azm/(r360/(r8*irrr)) iaaa=mod(iaaa,8*irrr) if(iaaa<0) iaaa=iaaa+8*irrr iaaa=iaaa+1 iaaamax=max(iaaamax,iaaa) iaaamin=min(iaaamin,iaaa) error = erradar_inflate*thiserr errmax=max(error,errmax) if(thiserr>zero) errmin=min(error,errmin) ! Perform limited qc based on azimuth angle, elevation angle, radial wind ! speed, range, distance from radar site good0=.true. if(abs(azm)>r400) then ibadazm=ibadazm+1; good0=.false. end if if(abs(rwnd) > r71) then ibadwnd=ibadwnd+1; good0=.false. end if if(thisrange>r92) then ibadrange=ibadrange+1; good0=.false. end if if(dist>r400) then ibaddist=ibaddist+1; good0=.false. end if if(height<-r1000.or.height>r50000) then ibadheight=ibadheight+1; good0=.false. end if good=.true. if(.not.good0) then notgood0=notgood0+1 cycle end if ! if data is good, load into output array if(good) then ntdrvr_kept=ntdrvr_kept+1 !#################### Data thinning ################### icntpnt=icntpnt+1 if(ithin > 0)then if(zflag == 0)then klon1= int(dlon); klat1= int(dlat) dx = dlon-klon1; dy = dlat-klat1 dx1 = one-dx; dy1 = one-dy w00=dx1*dy1; w10=dx1*dy; w01=dx*dy1; w11=dx*dy klat1=min(max(1,klat1),nlat); klon1=min(max(0,klon1),nlon) if (klon1==0) klon1=nlon klatp1=min(nlat,klat1+1); klonp1=klon1+1 if (klonp1==nlon+1) klonp1=1 do kk=1,nsig hges(kk)=w00*hgtl_full(klat1 ,klon1 ,kk) + & w10*hgtl_full(klatp1,klon1 ,kk) + & w01*hgtl_full(klat1 ,klonp1,kk) + & w11*hgtl_full(klatp1,klonp1,kk) end do sin2 = sin(dlat_earth)*sin(dlat_earth) termg = grav_equator * & ((one+somigliana*sin2)/sqrt(one-eccentricity*eccentricity*sin2)) termr = semi_major_axis /(one + flattening + grav_ratio - & two*flattening*sin2) termrg = (termg/grav)*termr do k=1,nsig zges(k) = (termr*hges(k)) / (termrg-hges(k)) zl_thin(k)=zges(k) end do endif zobs = height ntmp=ndata ! counting moved to map3gridS if (thin4d) then timedif = zero else timedif=abs(t4dv-toff) endif crit1 = timedif/r6+half call map3grids(1,zflag,zl_thin,nlevz,dlat_earth,dlon_earth,& zobs,crit1,ndata,iout,icntpnt,iiout,luse,.false.,.false.) maxout=max(maxout,iout) maxdata=max(maxdata,ndata) if (.not. luse) then ntdrvr_thin2=ntdrvr_thin2+1 cycle endif if(iiout > 0) isort(iiout)=0 if (ndata > ntmp) then nodata=nodata+1 endif isort(icntpnt)=iout else ndata =ndata+1 nodata=nodata+1 iout=ndata isort(icntpnt)=iout endif if(ndata > maxobs) then write(6,*)'READ_PREPBUFR: ***WARNING*** ndata > maxobs for ',obstype ndata = maxobs end if ! Set usage variable usage = zero if(icuse(ikx) < 0)usage=r100 if(ncnumgrp(ikx) > 0 )then ! cross validation on if(mod(ndata,ncnumgrp(ikx))== ncgroup(ikx)-1)usage=ncmiter(ikx) end if call deter_zsfc_model(dlat,dlon,zsges) ! Get information from surface file necessary for conventional data here call deter_sfc2(dlat_earth,dlon_earth,t4dv,idomsfc,skint,ff10,sfcr) cdata(1) = error ! wind obs error (m/s) cdata(2) = dlon ! grid relative longitude cdata(3) = dlat ! grid relative latitude cdata(4) = height ! obs absolute height (m) cdata(5) = rwnd ! wind obs (m/s) cdata(6) = azm*deg2rad ! azimuth angle (radians) cdata(7) = t4dv ! obs time (hour) cdata(8) = ikx ! type cdata(9) = tiltangle ! tilt angle (radians) cdata(10)= staheight ! station elevation (m) cdata(11)= rstation_id ! station id cdata(12)= usage ! usage parameter cdata(13)= idomsfc ! dominate surface type cdata(14)= skint ! skin temperature cdata(15)= ff10 ! 10 meter wind factor cdata(16)= sfcr ! surface roughness cdata(17)=dlon_earth_deg ! earth relative longitude (degrees) cdata(18)=dlat_earth_deg ! earth relative latitude (degrees) cdata(19)=dist ! range from radar in km (used to estimate beam spread) cdata(20)=zsges ! model elevation at radar site cdata(21)=thiserr cdata(22)=three+two ! tail Doppler radar do j=1,maxdat cdata_all(j,iout)=cdata(j) end do jjj=jjj+1 else notgood = notgood + 1 end if ! if(good) end do loop3! end of loop, reading records of data close(25) end do ! end of loop, reading TDR so data files close(lnbufr) else nswptype=0 nmrecs=0 irec=0 ! Open data file open(lnbufr,file=trim(infile),form='unformatted') call openbf(lnbufr,'IN',lnbufr) call datelen(10) call readmg(lnbufr,subset,idate,iret) if(iret==0) then ! Time offset call time_4dvar(idate,toff) write(date,'( i10)') idate read (date,'(i4,3i2)') iy,im,idd,ihh write(6,*)'READ_RADAR: bufr file date is ',iy,im,idd,ihh idate5(1) = iy ! year idate5(2) = im ! month idate5(3) = idd ! day idate5(4) = ihh ! hour idate5(5) = 0 ! minute call w3fs21(idate5,mincy) if(l_foreaft_thin)then ! Read the first 500 records to deterine which criterion ! should be used to seperate fore/aft sweep ! Big loop over bufr file loop5: do call readsb(lnbufr,iret) if(iret/=0) then call readmg(lnbufr,subset,idate,iret) if(iret/=0) exit loop5 cycle loop5 end if if(subset/=subset_check(loop)) then call readmg(lnbufr,subset,idate,iret) if(iret/=0) exit loop5 cycle loop5 end if nmrecs = nmrecs+1 ! Read header. Extract elevation angle call ufbint(lnbufr,hdr,12,1,levs,hdrstr(2)) thistilt=hdr(12) if(nmrecs == 1)then tdrele1 = hdr(12) tdrele2 = hdr(12) end if tdrele1 = tdrele2 tdrele2 = hdr(12) if(abs(tdrele2-tdrele1)>r100) then print *,'tdrele2,tdrele1=',tdrele2,tdrele1 nswptype=1 exit loop5 end if if(nmrecs <= 500)then cycle loop5 else exit loop5 end if end do loop5 firstbeam = 0 foreswp = .true. aftswp = .false. nforeswp=1 naftswp=0 nswp=1 call closbf(lnbufr) close(lnbufr) open(lnbufr,file=trim(infile),form='unformatted') call openbf(lnbufr,'IN',lnbufr) call datelen(10) call readmg(lnbufr,subset,idate,iret) else foreswp = .false. aftswp = .false. end if print *,'nmrecs, nswptype=', nmrecs, nswptype nmrecs=0 ! Big loop over bufr file loop4: do call readsb(lnbufr,iret) if(iret/=0) then call readmg(lnbufr,subset,idate,iret) if(iret/=0) exit loop4 cycle loop4 end if if(subset/=subset_check(loop)) then call readmg(lnbufr,subset,idate,iret) if(iret/=0) exit loop4 cycle loop4 end if nmrecs = nmrecs+1 irec = irec+1 ! Read header. Extract station infomration call ufbint(lnbufr,hdr,12,1,levs,hdrstr(2)) ! rstation_id=hdr(1) if(hdr(1) == zero)then cstaid='NOAA ' else if(hdr(1) == one)then cstaid='FRENCH ' else if(hdr(1)== two)then cstaid='G-IV ' else if(hdr(1)== three)then cstaid='AOC ' else cstaid='UNKNOWN ' endif kx=990+nint(hdr(1)) if(nmrecs==1)print *,'Antenna ID:', hdr(1),cstaid iyr = hdr(2) imo = hdr(3) idy = hdr(4) ihr = hdr(5) imn = hdr(6) isc = hdr(7) idate5(1) = iyr idate5(2) = imo idate5(3) = idy idate5(4) = ihr idate5(5) = imn ikx=0 do i=1,nconvtype if(trim(ioctype(i)) == trim(obstype) .and. kx == ictype(i))ikx = i end do if(ikx == 0) cycle loop4 call w3fs21(idate5,minobs) t4dv=real(minobs-iwinbgn,r_kind)*r60inv if (l4dvar.or.l4densvar) then if (t4dvwinlen) then ntimeout=ntimeout+1 cycle loop4 end if timeo=t4dv else timeo = real(minobs-mincy,r_kind)*r60inv if (abs(timeo) > twind .or. abs(timeo) > ctwind(ikx)) then ntimeout=ntimeout+1 cycle loop4 end if endif timemax=max(timemax,timeo) timemin=min(timemin,timeo) this_stalat=hdr(8) this_stalon=hdr(9) rlon0=deg2rad*this_stalon this_stalatr=this_stalat*deg2rad clat0=cos(this_stalatr) ; slat0=sin(this_stalatr) this_stahgt=hdr(10) thisazimuth=hdr(11) thistilt=hdr(12) elevmax=max(elevmax,thistilt) elevmin=min(elevmin,thistilt) ! define fore/aft sweeps for thinning (pseduo dual Doppler) if(l_foreaft_thin)then if (firstbeam == 0) then tdrele1 = hdr(12) tdrele2 = hdr(12) if(nswptype == 0)then tdrele3 = hdr(12) end if firstbeam = 1 endif if(nswptype == 0)then tdrele1 = tdrele2 tdrele2 = tdrele3 tdrele3 = hdr(12) if(firstbeam > 0 .and. tdrele2>=tdrele1 .and. tdrele2>=tdrele3 .and. tdrele2 > r60 & .and. irec > r150)then if(foreswp) then foreswp = .false. aftswp = .true. naftswp = naftswp+1 irec=0 else aftswp = .false. foreswp = .true. nforeswp = nforeswp+1 irec=0 endif nswp = nswp+1 endif else if(nswptype == 1)then tdrele1 = tdrele2 tdrele2 = hdr(12) if(abs(tdrele2-tdrele1)>r100) then if(foreswp) then foreswp = .false. aftswp = .true. naftswp = naftswp+1 irec=0 else aftswp = .false. foreswp = .true. nforeswp = nforeswp+1 irec=0 endif nswp = nswp+1 endif else foreswp = .false. aftswp = .false. end if else foreswp = .false. aftswp = .false. endif if(abs(thistilt)>r75)then ibadtilt=ibadtilt+1; cycle loop4 endif staheight=this_stahgt if(staheight<-r1000.or.staheight>r50000) then ibadstaheight=ibadstaheight+1; cycle loop4 end if ! Go through the data levels call ufbint(lnbufr,tdr_obs,4,maxlevs,levs,datstr(2)) if(levs>maxlevs) then write(6,*)'READ_RADAR: ***ERROR*** increase read_radar bufr size since ',& 'number of levs=',levs,' > maxlevs=',maxlevs call stop2(84) endif ! use local coordinate centered on this_stalat,this_stalon. note that global and local ! azimuth angle are the same at the origin (this_stalat,this_stalon) ! and azimuth angle is fixed in local coordinate along entire radial line. ! we convert back to global azimuth angle at each point along line ! at end of computation. that way we avoid worrying about where poles are. aactual=erad+this_stahgt thistiltr=thistilt*deg2rad selev0=sin(thistiltr) ; celev0=cos(thistiltr) a43=four_thirds*aactual ii=0 do k=1,levs nread=nread+1 ! Select data every 3 km along each beam if(MOD(INT(tdr_obs(1,k)-tdr_obs(1,1)),3000) < 100)then if(tdr_obs(3,k) >= 800.) nmissing=nmissing+1 !xx if(tdr_obs(3,k) < 800.) then ii=ii+1 dopbin(ii)=tdr_obs(3,k) thisrange=tdr_obs(1,k) call getvrlocalinfo(thisrange,thisazimuth,this_stahgt,aactual,a43,selev0,celev0, & rlon0,clat0,slat0,r8,r89_5,nsubzero,ii,z(ii),elev(ii),elat8(ii), & elon8(ii),glob_azimuth8(ii)) end if else ntdrvr_thin1=ntdrvr_thin1+1 endif end do ! Further process tail Doppler radar Vr data iimax=max(iimax,ii) if( ii > 0 )then dlat_earth=this_stalat !station lat (degrees) dlon_earth=this_stalon !station lon (degrees) if (dlon_earth>=r360) dlon_earth=dlon_earth-r360 if (dlon_earth=r360) dlon_earth=dlon_earth-r360 if(dlon_earthmax_rrr)then nirrr=nirrr+1 cycle endif ! Extract radial wind data height= z(i) rwnd = dopbin(i) azm_earth = glob_azimuth8(i) if(regional .and. .not. fv3_regional) then cosazm_earth=cos(azm_earth*deg2rad) sinazm_earth=sin(azm_earth*deg2rad) call rotate_wind_ll2xy(cosazm_earth,sinazm_earth,cosazm,sinazm,dlon_earth,dlon,dlat) azm=atan2(sinazm,cosazm)*rad2deg else azm=azm_earth end if iaaa=azm/(r360/(r8*irrr)) iaaa=mod(iaaa,8*irrr) if(iaaa<0) iaaa=iaaa+8*irrr iaaa=iaaa+1 iaaamax=max(iaaamax,iaaa) iaaamin=min(iaaamin,iaaa) error = erradar_inflate*thiserr errmax=max(error,errmax) if(thiserr>zero) errmin=min(error,errmin) ! Perform limited qc based on azimuth angle, elevation angle, radial wind ! speed, range, distance from radar site good0=.true. if(abs(azm)>r400) then ibadazm=ibadazm+1; good0=.false. end if if(abs(rwnd) > r71 .or. abs(rwnd) < r2 ) then ibadwnd=ibadwnd+1; good0=.false. end if if(thisrange>r92) then ibadrange=ibadrange+1; good0=.false. end if if(dist>r400) then ibaddist=ibaddist+1; good0=.false. end if if(height<-r1000.or.height>r50000) then ibadheight=ibadheight+1; good0=.false. end if good=.true. if(.not.good0) then notgood0=notgood0+1 cycle end if ! if data is good, load into output array if(good) then ntdrvr_kept=ntdrvr_kept+1 !#################### Data thinning ################### icntpnt=icntpnt+1 if(ithin > 0)then if(zflag == 0)then klon1= int(dlon); klat1= int(dlat) dx = dlon-klon1; dy = dlat-klat1 dx1 = one-dx; dy1 = one-dy w00=dx1*dy1; w10=dx1*dy; w01=dx*dy1; w11=dx*dy klat1=min(max(1,klat1),nlat); klon1=min(max(0,klon1),nlon) if (klon1==0) klon1=nlon klatp1=min(nlat,klat1+1); klonp1=klon1+1 if (klonp1==nlon+1) klonp1=1 do kk=1,nsig hges(kk)=w00*hgtl_full(klat1 ,klon1 ,kk) + & w10*hgtl_full(klatp1,klon1 ,kk) + & w01*hgtl_full(klat1 ,klonp1,kk) + & w11*hgtl_full(klatp1,klonp1,kk) end do sin2 = sin(dlat_earth)*sin(dlat_earth) termg = grav_equator * & ((one+somigliana*sin2)/sqrt(one-eccentricity*eccentricity*sin2)) termr = semi_major_axis /(one + flattening + grav_ratio - & two*flattening*sin2) termrg = (termg/grav)*termr do k=1,nsig zges(k) = (termr*hges(k)) / (termrg-hges(k)) zl_thin(k)=zges(k) end do endif zobs = height ntmp=ndata ! counting moved to map3gridS if (thin4d) then timedif = zero else timedif=abs(t4dv-toff) endif crit1 = timedif/r6+half call map3grids(1,zflag,zl_thin,nlevz,dlat_earth,dlon_earth,& zobs,crit1,ndata,iout,icntpnt,iiout,luse,foreswp,aftswp) maxout=max(maxout,iout) maxdata=max(maxdata,ndata) if (.not. luse) then if (foreswp) then ntdrvr_thin2_foreswp=ntdrvr_thin2_foreswp+1 else if (aftswp) then ntdrvr_thin2_aftswp=ntdrvr_thin2_aftswp+1 end if ntdrvr_thin2=ntdrvr_thin2+1 cycle endif if(iiout > 0) isort(iiout)=0 if (ndata > ntmp) then nodata=nodata+1 endif isort(icntpnt)=iout else ndata =ndata+1 nodata=nodata+1 iout=ndata isort(icntpnt)=iout endif if(ndata > maxobs) then write(6,*)'READ_PREPBUFR: ***WARNING*** ndata > maxobs for ',obstype ndata = maxobs end if ! Set usage variable usage = zero if(icuse(ikx) < 0)usage=r100 if(ncnumgrp(ikx) > 0 )then ! cross validation on if(mod(ndata,ncnumgrp(ikx))== ncgroup(ikx)-1)usage=ncmiter(ikx) end if call deter_zsfc_model(dlat,dlon,zsges) ! Get information from surface file necessary for conventional data here call deter_sfc2(dlat_earth,dlon_earth,t4dv,idomsfc,skint,ff10,sfcr) cdata(1) = error ! wind obs error (m/s) cdata(2) = dlon ! grid relative longitude cdata(3) = dlat ! grid relative latitude cdata(4) = height ! obs absolute height (m) cdata(5) = rwnd ! wind obs (m/s) cdata(6) = azm*deg2rad ! azimuth angle (radians) cdata(7) = t4dv ! obs time (hour) cdata(8) = ikx ! type cdata(9) = tiltangle ! tilt angle (radians) cdata(10)= staheight ! station elevation (m) cdata(11)= rstation_id ! station id cdata(12)= usage ! usage parameter cdata(13)= idomsfc ! dominate surface type cdata(14)= skint ! skin temperature cdata(15)= ff10 ! 10 meter wind factor cdata(16)= sfcr ! surface roughness cdata(17)=dlon_earth_deg ! earth relative longitude (degrees) cdata(18)=dlat_earth_deg ! earth relative latitude (degrees) cdata(19)=dist ! range from radar in km (used to estimate beam spread) cdata(20)=zsges ! model elevation at radar site cdata(21)=thiserr cdata(22)=hdr(1)+three+one ! tail Doppler radar do j=1,maxdat cdata_all(j,iout)=cdata(j) end do if(foreswp)nfore=nfore+1 if(aftswp)naft=naft+1 jjj=jjj+1 else notgood = notgood + 1 end if ! if(good) end do endif ! if(ii .gt. 0) ! End of bufr read loop end do loop4 ! Normal exit else write(6,*)'READ_RADAR: problem reading tail Doppler radar bufr file tldplrbufr' end if call closbf(lnbufr) end if 1200 continue close(lnbufr) if (.not. use_all) then deallocate(zl_thin) call del3grids endif write(6,*)'READ_RADAR: # records(beams) read in nmrecs=', nmrecs write(6,*)'READ_RADAR: # records out of time window =', ntimeout write(6,*)'READ_RADAR: # records with bad tilt=',ibadtilt write(6,*)'READ_RADAR: # records with bad station height =',ibadstaheight write(6,*)'READ_RADAR: # data read in nread=', nread write(6,*)'READ_RADAR: # data with missing value nmissing=', nmissing write(6,*)'READ_RADAR: # data likely to be below sealevel nsubzero=', nsubzero write(6,*)'READ_RADAR: # data removed by thinning along the beam ntdrvr_thin1=', ntdrvr_thin1 write(6,*)'READ_RADAR: # data retained after thinning along the beam ntdrvr_in=', ntdrvr_in write(6,*)'READ_RADAR: # out of domain =', noutside write(6,*)'READ_RADAR: # out of range =', nirrr write(6,*)'READ_RADAR: # bad azimuths =',ibadazm write(6,*)'READ_RADAR: # bad winds (<2m/s or >71m/s) =',ibadwnd write(6,*)'READ_RADAR: # bad ranges =',ibadrange write(6,*)'READ_RADAR: # bad distance from radar =',ibaddist write(6,*)'READ_RADAR: # bad obs height =',ibadheight write(6,*)'READ_RADAR: # bad data =',notgood0 write(6,*)'READ_RADAR: # data retained after QC ntdrvr_kept=', ntdrvr_kept write(6,*)'READ_RADAR: # data removed by thinning mesh ntdrvr_thin2=', ntdrvr_thin2 if(l_foreaft_thin)then write(6,*)'READ_RADAR: nforeswp,naftswp,nswp=',nforeswp,naftswp,nswp write(6,*)'READ_RADAR: ntdrvr_thin2_foreswp,ntdrvr_thin2_aftswp=',ntdrvr_thin2_foreswp,ntdrvr_thin2_aftswp write(6,*)'READ_RADAR: data retained for further processing nfore,naft=',nfore,naft end if write(6,*)'READ_RADAR: data retained for further processing =', jjj write(6,*)'READ_RADAR: timemin,max =',timemin,timemax write(6,*)'READ_RADAR: elevmin,max =',elevmin,elevmax write(6,*)'READ_RADAR: dlatmin,max,dlonmin,max=',dlatmin,dlatmax,dlonmin,dlonmax write(6,*)'READ_RADAR: iaaamin,max,8*max_rrr =',iaaamin,iaaamax,8*max_rrr write(6,*)'READ_RADAR: iimax =',iimax ! Write observation to scratch file call count_obs(ndata,maxdat,ilat,ilon,cdata_all,nobs) write(lunout) obstype,sis,nreal,nchanl,ilat,ilon write(lunout) ((cdata_all(k,i),k=1,maxdat),i=1,ndata) deallocate(cdata_all) return 300 write(6,*) 'read_radar open TDR SO file list failed ' call stop2(555) 200 write(6,*) 'read_radar read TDR SO data failed ' call stop2(555) end subroutine read_radar subroutine getvrlocalinfo(thisrange,thisazimuth,this_stahgt,aactual,a43,selev0,celev0, & rlon0,clat0,slat0,r8,r89_5,nsubzero,ii,z,elev,elat8,elon8, & glob_azimuth8) !$$$ subprogram documentation block ! . . . . ! subprogram: getvrlocalinfo following subroutine radar_bufr_read_all ! prgmmr: tong org: np23 date: 2013-03-28 ! ! abstract: This routine calcuate radial wind elevation, elevation angle, ! earth lat lon and and azimuth angle at observation location ! ! program history log: ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block use kinds, only: r_kind,r_single,i_kind use constants, only: one,half,two,deg2rad,rad2deg,zero_single,rearth use read_l2bufr_mod, only: invtllv implicit none real(r_single) ,intent(in ) :: thisrange,thisazimuth,a43,aactual,selev0,celev0 real(r_kind) ,intent(in ) :: this_stahgt,rlon0,clat0,slat0,r8,r89_5 integer(i_kind),intent(inout) :: nsubzero integer(i_kind),intent(inout) :: ii real(r_single) ,intent(out ) :: elev,z,elat8,elon8,glob_azimuth8 ! local variables real(r_single) b,c,epsh,h,ha,celev,selev,gamma real(r_single) rad_per_meter real(r_kind) thisazimuthr,rlonloc,rlatloc,rlonglob,rlatglob,thislat,thislon real(r_kind) clat1,caz0,saz0,cdlon,sdlon,caz1,saz1 rad_per_meter= one/rearth ! use 4/3rds rule to get elevation of radar beam ! (if local temperature available, then vertical position can be ! estimated with greater accuracy) b=thisrange*(thisrange+two*aactual*selev0) c=sqrt(aactual*aactual+b) ha=b/(aactual+c) epsh=(thisrange*thisrange-ha*ha)/(r8*aactual) h=ha-epsh z=this_stahgt+h if(z < zero_single)then ! don't use observation if it is likely to be below sealevel nsubzero=nsubzero+1 ii=ii-1 else ! Get elevation angle at obs location celev=celev0 selev=selev0 if(thisrange>=one) then celev=a43*celev0/(a43+h) selev=(thisrange*thisrange+h*h+two*a43*h)/(two*thisrange*(a43+h)) end if elev=rad2deg*atan2(selev,celev) gamma=half*thisrange*(celev0+celev) ! Get earth lat lon at obs location thisazimuthr=thisazimuth*deg2rad rlonloc=rad_per_meter*gamma*cos(thisazimuthr) rlatloc=rad_per_meter*gamma*sin(thisazimuthr) call invtllv(rlonloc,rlatloc,rlon0,clat0,slat0,rlonglob,rlatglob) thislat=rlatglob*rad2deg thislon=rlonglob*rad2deg ! Keep away from poles if(abs(thislat)>r89_5)then ii=ii-1 else elat8=thislat elon8=thislon ! Get corrected azimuth clat1=cos(rlatglob) caz0=cos(thisazimuthr) saz0=sin(thisazimuthr) cdlon=cos(rlonglob-rlon0) sdlon=sin(rlonglob-rlon0) caz1=clat0*caz0/clat1 saz1=saz0*cdlon-caz0*sdlon*slat0 glob_azimuth8=atan2(saz1,caz1)*rad2deg end if end if return end subroutine getvrlocalinfo subroutine read_radar_l2rw_novadqc(ndata,nodata,lunout,obstype,sis,nobs) !$$$ subprogram documentation block ! . . . . ! subprogram: read_radar_l2rw_novadqc read radar L2 radial winds no VAD QC ! prgmmr: yang org: np23 date: 1998-05-15 ! ! abstract: This routine reads radar radial wind files. ! ! When running the gsi in regional mode, the code only ! retains those observations that fall within the regional ! ! program history log: ! 2015-10-19 lippi - Modified from read_radar to only process level 2 radial ! wind obs. and skip vad wind checks. ! ! input argument list: ! lunout - unit to which to write data for further processing ! obstype - observation type to process ! ! output argument list: ! ndata - number of doppler lidar wind profiles retained for further ! processing ! nodata - number of doppler lidar wind observations retained for further ! processing ! sis - satellite/instrument/sensor indicator ! nobs - array of observations on each subdomain for each processor! use kinds, only: r_kind,r_single,r_double,i_kind,i_byte use constants, only: zero,half,one,two,deg2rad,rearth,rad2deg,r1000,r100,r400 use qcmod, only: erradar_inflate use oneobmod, only: oneobtest,learthrel_rw use gsi_4dvar, only: l4dvar,l4densvar,winlen,time_4dvar use gridmod, only: regional,nlat,nlon,tll2xy,rlats,rlons,rotate_wind_ll2xy,& fv3_regional use convinfo, only: nconvtype,ncmiter,ncgroup,ncnumgrp,icuse,ioctype use deter_sfc_mod, only: deter_sfc2 use mpimod, only: npe implicit none ! Declare passed variables character(len=*),intent(in ) :: obstype!,infile character(len=20),intent(in ) :: sis ! real(r_kind) ,intent(in ) :: twind integer(i_kind) ,intent(in ) :: lunout integer(i_kind) ,intent(inout) :: ndata,nodata!,nread integer(i_kind),dimension(npe) ,intent(inout) :: nobs ! Declare local parameters integer(i_kind),parameter:: maxlevs=1500 integer(i_kind),parameter:: maxdat=22 real(r_kind),parameter:: r4_r_kind = 4.0_r_kind real(r_kind),parameter:: r6 = 6.0_r_kind real(r_kind),parameter:: r8 = 8.0_r_kind real(r_kind),parameter:: r90 = 90.0_r_kind real(r_kind),parameter:: r200 = 200.0_r_kind real(r_kind),parameter:: r150 = 150.0_r_kind real(r_kind),parameter:: r360 = 360.0_r_kind real(r_kind),parameter:: r50000 = 50000.0_r_kind real(r_kind),parameter:: r89_5 = 89.5_r_kind real(r_kind),parameter:: four_thirds = 4.0_r_kind / 3.0_r_kind ! Declare local variables logical good,outside,good0 character(30) outmessage integer(i_kind) lnbufr,i,k,maxobs integer(i_kind) nmrecs,ibadazm,ibadwnd,ibaddist,ibadheight,kthin integer(i_kind) ibadstaheight,ibaderror,notgood,iheightbelowsta,ibadfit integer(i_kind) notgood0 integer(i_kind) iret,kx0 integer(i_kind) nreal,nchanl,ilat,ilon,ikx integer(i_kind) idomsfc real(r_kind) usage,ff10,sfcr,skint,t4dv,t4dvo,toff real(r_kind) eradkm,dlat_earth,dlon_earth real(r_kind) dlat,dlon,staheight,tiltangle,clon,slon,clat,slat real(r_kind) timeo,clonh,slonh,clath,slath,cdist,dist real(r_kind) rwnd,azm,height,error real(r_kind) azm_earth,cosazm_earth,sinazm_earth,cosazm,sinazm real(r_kind):: zsges real(r_kind),dimension(maxdat):: cdata real(r_kind),allocatable,dimension(:,:):: cdata_all real(r_double) rstation_id character(8) cstaid character(4) this_staid equivalence (this_staid,cstaid) equivalence (cstaid,rstation_id) integer(i_kind) loop real(r_kind) timemax,timemin,errmax,errmin real(r_kind) dlatmax,dlonmax,dlatmin,dlonmin real(r_kind) xscale,xscalei integer(i_kind) max_rrr,nboxmax integer(i_kind) irrr,iaaa,iaaamax,iaaamin real(r_kind) this_stalat,this_stalon,this_stahgt,thistime,thislat,thislon real(r_kind) thishgt,thisvr,corrected_azimuth,thiserr,corrected_tilt integer(i_kind) nsuper2_in,nsuper2_kept real(r_kind) errzmax integer(i_kind),allocatable,dimension(:):: isort ! following variables are for fore/aft separation integer(i_kind) irec data lnbufr/10/ !*********************************************************************************** eradkm=rearth*0.001_r_kind maxobs=2e6 nreal=maxdat nchanl=0 ilon=2 ilat=3 iaaamax=-huge(iaaamax) iaaamin=huge(iaaamin) dlatmax=-huge(dlatmax) dlonmax=-huge(dlonmax) dlatmin=huge(dlatmin) dlonmin=huge(dlonmin) allocate(cdata_all(maxdat,maxobs),isort(maxobs)) isort = 0 cdata_all=zero ! Initialize variables xscale=1000._r_kind xscalei=one/xscale max_rrr=nint(200000.0_r_kind*xscalei) nboxmax=1 kx0=22500 nmrecs=0 irec=0 errzmax=zero ! First process any level 2 superobs. ! Initialize variables. ikx=0 do i=1,nconvtype if(trim(ioctype(i)) == trim(obstype))ikx = i end do timemax=-huge(timemax) timemin=huge(timemin) errmax=-huge(errmax) errmin=huge(errmin) loop=0 ibadazm=0 ibadwnd=0 ibaddist=0 ibadheight=0 ibadstaheight=0 iheightbelowsta=0 iheightbelowsta=0 ibaderror=0 ibadfit=0 kthin=0 notgood=0 notgood0=0 nsuper2_in=0 nsuper2_kept=0 if(loop==0) outmessage='level 2 superobs:' ! Open sequential file containing superobs open(lnbufr,file='radar_supobs_from_level2',form='unformatted') rewind lnbufr ! Loop to read superobs data file do read(lnbufr,iostat=iret)this_staid,this_stalat,this_stalon,this_stahgt, & thistime,thislat,thislon,thishgt,thisvr,corrected_azimuth,thiserr,corrected_tilt if(iret/=0) exit nsuper2_in=nsuper2_in+1 dlat_earth=this_stalat !station lat (degrees) dlon_earth=this_stalon !station lon (degrees) if (dlon_earth>=r360) dlon_earth=dlon_earth-r360 if (dlon_earthwinlen) cycle else timeo=thistime if(abs(timeo)>half ) cycle endif ! Get observation (lon,lat). Compute distance from radar. dlat_earth=thislat dlon_earth=thislon if(dlon_earth>=r360) dlon_earth=dlon_earth-r360 if(dlon_earthmax_rrr) cycle end if ! Extract radial wind data height= thishgt rwnd = thisvr azm_earth = corrected_azimuth if(regional .and. .not. fv3_regional) then if(oneobtest .and. learthrel_rw) then ! for non rotated winds!!! cosazm=cos(azm_earth*deg2rad) sinazm=sin(azm_earth*deg2rad) azm=atan2(sinazm,cosazm)*rad2deg else cosazm_earth=cos(azm_earth*deg2rad) sinazm_earth=sin(azm_earth*deg2rad) call rotate_wind_ll2xy(cosazm_earth,sinazm_earth,cosazm,sinazm,dlon_earth,dlon,dlat) azm=atan2(sinazm,cosazm)*rad2deg end if else azm=azm_earth end if if(.not. oneobtest) then iaaa=azm/(r360/(r8*irrr)) iaaa=mod(iaaa,8*irrr) if(iaaa<0) iaaa=iaaa+8*irrr iaaa=iaaa+1 iaaamax=max(iaaamax,iaaa) iaaamin=min(iaaamin,iaaa) end if error = erradar_inflate*thiserr errmax=max(error,errmax) if(thiserr>zero) errmin=min(error,errmin) ! Perform limited qc based on azimuth angle, radial wind ! speed, distance from radar site, elevation of radar, ! height of observation, and observation error good0=.true. if(abs(azm)>r400) then ibadazm=ibadazm+1; good0=.false. end if if(abs(rwnd)>r200) then ibadwnd=ibadwnd+1; good0=.false. end if if(dist>r400) then ibaddist=ibaddist+1; good0=.false. end if if(staheight<-r1000.or.staheight>r50000) then ibadstaheight=ibadstaheight+1; good0=.false. end if if(height<-r1000.or.height>r50000) then ibadheight=ibadheight+1; good0=.false. end if if(heightr6 .or. thiserr<=zero) then ibaderror=ibaderror+1; good0=.false. end if good=.true. if(.not.good0) then notgood0=notgood0+1 cycle else end if ! If data is good, load into output array if(good) then nsuper2_kept=nsuper2_kept+1 ndata =min(ndata+1,maxobs) nodata =min(nodata+1,maxobs) !number of obs not used (no meaninghere) usage = zero if(icuse(ikx) < 0)usage=r100 if(ncnumgrp(ikx) > 0 )then ! cross validation on if(mod(ndata,ncnumgrp(ikx))== ncgroup(ikx)-1)usage=ncmiter(ikx) end if call deter_sfc2(dlat_earth,dlon_earth,t4dv,idomsfc,skint,ff10,sfcr) cdata(1) = error ! wind obs error (m/s) cdata(2) = dlon ! grid relative longitude cdata(3) = dlat ! grid relative latitude cdata(4) = height ! obs absolute height (m) cdata(5) = rwnd ! wind obs (m/s) cdata(6) = azm*deg2rad ! azimuth angle (radians) cdata(7) = t4dv ! obs time (hour) cdata(8) = ikx ! type cdata(9) = tiltangle ! tilt angle (radians) cdata(10)= staheight ! station elevation (m) cdata(11)= rstation_id ! station id cdata(12)= usage ! usage parameter cdata(13)= idomsfc ! dominate surface type cdata(14)= skint ! skin temperature cdata(15)= ff10 ! 10 meter wind factor cdata(16)= sfcr ! surface roughness cdata(17)=dlon_earth*rad2deg ! earth relative longitude (degrees) cdata(18)=dlat_earth*rad2deg ! earth relative latitude (degrees) cdata(19)=dist ! range from radar in km (used to estimatebeam spread) cdata(20)=zsges ! model elevation at radar site cdata(21)=thiserr cdata(22)=two do i=1,maxdat cdata_all(i,ndata)=cdata(i) end do else notgood = notgood + 1 end if end do close(lnbufr) ! A simple unformatted fortran file should not be mixed with bufr I/O write(6,*)'READ_RADAR_L2RW_NOVADQC: ',trim(outmessage),' reached eof on 2 superob radar file' write(6,*)'READ_RADAR_L2RW_NOVADQC: nsuper2_in,nsuper2_kept=',nsuper2_in,nsuper2_kept write(6,*)'READ_RADAR_L2RW_NOVADQC: # bad azimuths=',ibadazm write(6,*)'READ_RADAR_L2RW_NOVADQC: # bad winds =',ibadwnd write(6,*)'READ_RADAR_L2RW_NOVADQC: # bad dists =',ibaddist write(6,*)'READ_RADAR_L2RW_NOVADQC: # bad stahgts =',ibadstaheight write(6,*)'READ_RADAR_L2RW_NOVADQC: # bad obshgts =',ibadheight write(6,*)'READ_RADAR_L2RW_NOVADQC: # bad errors =',ibaderror write(6,*)'READ_RADAR_L2RW_NOVADQC: # bad fit =',ibadfit write(6,*)'READ_RADAR_L2RW_NOVADQC: # num thinned =',kthin write(6,*)'READ_RADAR_L2RW_NOVADQC: # notgood0 =',notgood0 write(6,*)'READ_RADAR_L2RW_NOVADQC: # notgood =',notgood write(6,*)'READ_RADAR_L2RW_NOVADQC: # hgt belowsta=',iheightbelowsta write(6,*)'READ_RADAR_L2RW_NOVADQC: timemin,max =',timemin,timemax write(6,*)'READ_RADAR_L2RW_NOVADQC: errmin,max =',errmin,errmax write(6,*)'READ_RADAR_L2RW_NOVADQC: dlatmin,max,dlonmin,max=',dlatmin,dlatmax,dlonmin,dlonmax write(6,*)'READ_RADAR_L2RW_NOVADQC: iaaamin,max,8*max_rrr=',iaaamin,iaaamax,8*max_rrr ! Write observation to scratch file call count_obs(ndata,maxdat,ilat,ilon,cdata_all,nobs) write(lunout) obstype,sis,nreal,nchanl,ilat,ilon write(lunout) ((cdata_all(k,i),k=1,maxdat),i=1,ndata) deallocate(cdata_all) return end subroutine read_radar_l2rw_novadqc !!!!!!!!!!!!!!! Added for l2rw thinning !!!!!!!!!!!!!!! subroutine read_radar_l2rw(ndata,nodata,lunout,obstype,sis,nobs,hgtl_full) use kinds, only: r_kind,r_single,r_double,i_kind,i_byte use constants, only: zero,half,one,two,deg2rad,rearth,rad2deg,r1000,r100,r400 use qcmod, only: erradar_inflate use oneobmod, only: oneobtest,learthrel_rw use gsi_4dvar, only: l4dvar,l4densvar,winlen,time_4dvar use gridmod, only: regional,nlat,nlon,tll2xy,rlats,rlons,rotate_wind_ll2xy,nsig use obsmod, only: doradaroneob,oneobradid,time_offset use mpeu_util, only: gettablesize,gettable use convinfo, only: nconvtype,icuse,ioctype use deter_sfc_mod, only: deter_sfc2 use mpimod, only: npe use read_l2bufr_mod, only: radar_sites,radar_rmesh,radar_zmesh,elev_angle_max,del_time,range_max use constants, only: eccentricity,somigliana,grav_ratio,grav,semi_major_axis,flattening,grav_equator use obsmod,only: radar_no_thinning,iadate use convthin, only: make3grids,map3grids implicit none ! Declare passed variables character(len=*),intent(in ) :: obstype character(len=20),intent(in ) :: sis integer(i_kind) ,intent(in ) :: lunout integer(i_kind) ,intent(inout) :: ndata,nodata integer(i_kind),dimension(npe) ,intent(inout) :: nobs real(r_kind),dimension(nlat,nlon,nsig),intent(in):: hgtl_full ! Declare local parameters integer(i_kind),parameter:: maxlevs=1500 integer(i_kind),parameter:: maxdat=22 real(r_kind),parameter:: r4_r_kind = 4.0_r_kind real(r_kind),parameter:: r6 = 6.0_r_kind real(r_kind),parameter:: r8 = 8.0_r_kind real(r_kind),parameter:: r90 = 90.0_r_kind real(r_kind),parameter:: r200 = 200.0_r_kind real(r_kind),parameter:: r150 = 150.0_r_kind real(r_kind),parameter:: r360 = 360.0_r_kind real(r_kind),parameter:: r50000 = 50000.0_r_kind real(r_kind),parameter:: r89_5 = 89.5_r_kind real(r_kind),parameter:: four_thirds = 4.0_r_kind / 3.0_r_kind integer(i_kind),parameter:: n_gates_max=4000 real(r_double),parameter:: r1e5_double = 1.0e5_r_double real(r_kind),parameter:: rinv60 = 1.0_r_kind/60.0_r_kind logical good,outside,good0 character(30) outmessage integer(i_kind) lnbufr,i,k,maxobs integer(i_kind) nmrecs,ibadazm,ibadwnd,ibaddist,ibadheight,kthin integer(i_kind) ibadstaheight,ibaderror,notgood,iheightbelowsta,ibadfit integer(i_kind) notgood0 integer(i_kind) iret,kx0 integer(i_kind) nreal,nchanl,ilat,ilon,ikx integer(i_kind) idomsfc real(r_kind) usage,ff10,sfcr,skint,t4dvo,toff real(r_kind) eradkm,dlat_earth,dlon_earth real(r_kind) dlat,dlon,staheight,tiltangle,clon,slon,clat,slat real(r_kind) timeo,clonh,slonh,clath,slath,cdist,dist real(r_kind) rwnd,azm,height,error real(r_kind) azm_earth,cosazm_earth,sinazm_earth,cosazm,sinazm real(r_kind):: zsges real(r_kind),dimension(maxdat):: cdata real(r_kind),allocatable,dimension(:,:):: cdata_all real(r_double) rstation_id character(8) cstaid character(4) this_staid equivalence (this_staid,cstaid) equivalence (cstaid,rstation_id) integer(i_kind) loop real(r_kind) timemax,timemin,errmax,errmin real(r_kind) dlatmax,dlonmax,dlatmin,dlonmin real(r_kind) xscale,xscalei integer(i_kind) max_rrr,nboxmax integer(i_kind) irrr,iaaa,iaaamax,iaaamin real(r_kind) this_stalat,this_stalon,this_stahgt,thistime,thislat,thislon real(r_kind) thishgt,thisvr,corrected_azimuth,thiserr,corrected_tilt integer(i_kind) nsuper2_in,nsuper2_kept real(r_kind) errzmax character(len=*),parameter:: tbname='SUPEROB_RADAR::' integer(i_kind) ntot,radar_true,radar_count,inbufr,lundx,idups,idate,n_gates,levs integer(i_kind) idate5(5) integer(i_kind) nminref,nminthis,nrange_max integer(i_kind) nobs_in,nradials_in,nradials_fail_angmax,nradials_fail_time,nradials_fail_elb,ireadmg,ireadsb integer(i_kind) nobs_badvr,nobs_badsr,j real(r_kind) rlon0,clat0,slat0,this_stalatr,thisrange,thisazimuth,thistilt,thisvr2 real(r_kind) rad_per_meter,erad,ddiffmin,distfact character(len=256),allocatable,dimension(:):: rtable character(4),allocatable,dimension(:):: rsite integer(i_kind),allocatable,dimension(:):: ruse character(8) chdr2,subset real(r_double) rdisttest(n_gates_max),hdr(10),hdr2(12),rwnd0(3,n_gates_max) character(4) stn_id equivalence (chdr2,hdr2(1)) real(r_kind) stn_lat,stn_lon,stn_hgt,stn_az,stn_el,t,range,vrmax,vrmin,aactual,a43,b,c,selev0,celev0,thistiltr,epsh,h,ha,rlonloc,rlatloc real(r_kind) celev,selev,gamma,thisazimuthr,rlonglob,rlatglob,clat1,caz0,saz0,cdlon,sdlon,caz1,saz1 real(r_kind):: relm,srlm,crlm,sph,cph,cc,anum,denom real(r_kind) :: rmesh,xmesh,zmesh,dx,dy,dx1,dy1,w00,w01,w10,w11 real(r_kind), allocatable, dimension(:) :: zl_thin integer(i_kind) :: ithin,zflag,nlevz,icntpnt,klon1,klat1,kk,klatp1,klonp1 real(r_kind),dimension(nsig):: hges,zges real(r_kind) sin2,termg,termr,termrg,zobs integer(i_kind) ntmp,iout,iiout,ntdrvr_thin2 real(r_kind) crit1,timedif integer(i_kind) maxout,maxdata logical :: luse integer(i_kind) iyref,imref,idref,ihref,nout integer(i_kind),allocatable,dimension(:):: isort ! following variables are for fore/aft separation integer(i_kind) irec data lnbufr/10/ if (radar_sites) then open(666,file=trim('gsiparm.anl'),form='formatted') call gettablesize(tbname,666,ntot,radar_count) allocate(rtable(radar_count),rsite(radar_count),ruse(radar_count)) call gettable(tbname,666,ntot,radar_count,rtable) do i=1,radar_count read(rtable(i),*) rsite(i),ruse(i) write(*,'(A14,X,A4,X,I3)'),"Radar Sites: ",rsite(i),ruse(i) end do end if rad_per_meter= one/rearth erad = rearth eradkm=rearth*0.001_r_kind maxobs=2e7 nreal=maxdat nchanl=0 ilon=2 ilat=3 ikx=0 do j=1,nconvtype if(trim(ioctype(j)) == trim(obstype))ikx = j end do iaaamax=-huge(iaaamax) iaaamin=huge(iaaamin) dlatmax=-huge(dlatmax) dlonmax=-huge(dlonmax) dlatmin=huge(dlatmin) dlonmin=huge(dlonmin) allocate(cdata_all(maxdat,maxobs),isort(maxobs)) isort = 0 cdata_all=zero xscale=1000._r_kind xscalei=one/xscale max_rrr=nint(1000000.0_r_kind*xscalei) nboxmax=1 kx0=22500 nmrecs=0 irec=0 errzmax=zero timemax=-huge(timemax) timemin=huge(timemin) errmax=-huge(errmax) errmin=huge(errmin) loop=0 ibadazm=0 ibadwnd=0 ibaddist=0 ibadheight=0 ibadstaheight=0 iheightbelowsta=0 iheightbelowsta=0 ibaderror=0 ibadfit=0 kthin=0 notgood=0 notgood0=0 nsuper2_in=0 nsuper2_kept=0 ntdrvr_thin2=0 maxout=0 maxdata=0 isort=0 icntpnt=0 nout=0 if(loop==0) outmessage='level 2 superobs:' rmesh=radar_rmesh zmesh=radar_zmesh nlevz=nint(16000._r_kind/zmesh) xmesh=rmesh call make3grids(xmesh,nlevz) allocate(zl_thin(nlevz)) zflag=1 if (zflag == 1) then do k=1,nlevz zl_thin(k)=k*zmesh enddo endif inbufr=10 open(inbufr,file="l2rwbufr",form='unformatted') rewind inbufr lundx=inbufr call openbf(inbufr,'IN',lundx) call datelen(10) iyref=iadate(1) imref=iadate(2) idref=iadate(3) ihref=iadate(4) idate5(1)=iyref idate5(2)=imref idate5(3)=idref idate5(4)=ihref idate5(5)=0 ! minutes call w3fs21(idate5,nminref) idups=0 nobs_in=0 nradials_in=0 nradials_fail_angmax=0 nradials_fail_time=0 nradials_fail_elb=0 ddiffmin=huge(ddiffmin) do while(ireadmg(inbufr,subset,idate)>=0) do while (ireadsb(inbufr)==0) call ufbint(inbufr,rdisttest,1,n_gates_max,n_gates,'DIST125M') if(n_gates>1) then do i=1,n_gates-1 if(nint(abs(rdisttest(i+1)-rdisttest(i)))==0) then idups=idups+1 else ddiffmin=min(abs(rdisttest(i+1)-rdisttest(i)),ddiffmin) end if end do end if distfact=zero if(nint(ddiffmin)==1) distfact=250._r_kind if(nint(ddiffmin)==2) distfact=125._r_kind if(distfact==zero) then write(6,*)'RADAR_BUFR_READ_ALL: problem with level 2 bufr file, gate distance scale factor undetermined, going with 125' distfact=125._r_kind end if call ufbint(inbufr,hdr2,12,1,levs,'SSTN CLAT CLON HSMSL HSALG ANEL YEAR MNTH DAYS HOUR MINU SECO') if(hdr2(6)>elev_angle_max) then nradials_fail_angmax=nradials_fail_angmax+1 cycle end if idate5(1)=nint(hdr2(7)) ; idate5(2)=nint(hdr2(8)) ; idate5(3)=nint(hdr2(9)) idate5(4)=nint(hdr2(10)) ; idate5(5)=nint(hdr2(11)) call w3fs21(idate5,nminthis) t=(real(nminthis-nminref,r_kind)+real(nint(hdr2(12)),r_kind)*rinv60)*rinv60 timemax=max(t,timemax) timemin=min(t,timemin) if(abs(t)>del_time) then nradials_fail_time=nradials_fail_time+1 cycle end if nobs_in=nobs_in+n_gates stn_id=chdr2 radar_true=0 if (radar_sites) then do i=1,radar_count if (trim(stn_id) .eq. trim(rsite(i)) .and. ruse(i) .eq. 1 ) radar_true=1 end do if (radar_true .eq. 0) cycle end if stn_lat=hdr2(2) stn_lon=hdr2(3) stn_hgt=hdr2(4)+hdr2(5) call ufbint(inbufr,hdr,10,1,levs, & 'SSTN YEAR MNTH DAYS HOUR MINU SECO ANAZ ANEL QCRW') nradials_in=nradials_in+1 stn_az=r90-hdr(8) stn_el=hdr(9) call ufbint(inbufr,rwnd0,3,n_gates_max,n_gates,'DIST125M DMVR DVSW') do i=1,n_gates range=distfact*rwnd0(1,i) if(range>range_max) then nrange_max=nrange_max+1 cycle end if if(rwnd0(2,i)>r1e5_double) then nobs_badvr=nobs_badvr+1 cycle end if if(rwnd0(3,i)>r1e5_double) then nobs_badsr=nobs_badsr+1 cycle end if this_stalat=stn_lat if(abs(this_stalat)>r89_5) cycle this_stalon=stn_lon rlon0=deg2rad*this_stalon this_stalatr=this_stalat*deg2rad clat0=cos(this_stalatr) ; slat0=sin(this_stalatr) this_staid=stn_id this_stahgt=stn_hgt thisrange= range thisazimuth=stn_az thistilt=stn_el thisvr=rwnd0(2,i) vrmax=max(vrmax,thisvr) vrmin=min(vrmin,thisvr) thisvr2=rwnd0(2,i)**2 thiserr=5.0_r_kind errmax=max(errmax,thiserr) errmin=min(errmin,thiserr) thistime=t aactual=erad+this_stahgt a43=four_thirds*aactual thistiltr=thistilt*deg2rad selev0=sin(thistiltr) celev0=cos(thistiltr) b=thisrange*(thisrange+two*aactual*selev0) c=sqrt(aactual*aactual+b) ha=b/(aactual+c) epsh=(thisrange*thisrange-ha*ha)/(r8*aactual) h=ha-epsh thishgt=this_stahgt+h celev=celev0 selev=selev0 if(thisrange>=one) then celev=a43*celev0/(a43+h) selev=(thisrange*thisrange+h*h+two*a43*h)/(two*thisrange*(a43+h)) end if corrected_tilt=atan2(selev,celev)*rad2deg gamma=half*thisrange*(celev0+celev) ! Get earth lat lon of superob thisazimuthr=thisazimuth*deg2rad rlonloc=rad_per_meter*gamma*cos(thisazimuthr) rlatloc=rad_per_meter*gamma*sin(thisazimuthr) RELM=rlonloc SRLM=SIN(RELM) CRLM=COS(RELM) SPH=SIN(rlatloc) CPH=COS(rlatloc) CC=CPH*CRLM ANUM=CPH*SRLM DENOM=clat0*CC-slat0*SPH rlonglob=rlon0+ATAN2(ANUM,DENOM) rlatglob=ASIN(clat0*SPH+slat0*CC) thislat=rlatglob*rad2deg thislon=rlonglob*rad2deg if(abs(thislat)>r89_5) cycle clat1=cos(rlatglob) caz0=cos(thisazimuthr) saz0=sin(thisazimuthr) cdlon=cos(rlonglob-rlon0) sdlon=sin(rlonglob-rlon0) caz1=clat0*caz0/clat1 saz1=saz0*cdlon-caz0*sdlon*slat0 corrected_azimuth=atan2(saz1,caz1)*rad2deg if (doradaroneob .and. (oneobradid /= this_staid)) cycle if(iret/=0) exit nsuper2_in=nsuper2_in+1 dlat_earth=this_stalat !station lat (degrees) dlon_earth=this_stalon !station lon (degrees) if (dlon_earth>=r360) dlon_earth=dlon_earth-r360 if (dlon_earthwinlen) cycle else timeo=thistime if(abs(timeo)>half ) cycle endif ! Get observation (lon,lat). Compute distance from radar. dlat_earth=thislat dlon_earth=thislon if(dlon_earth>=r360) dlon_earth=dlon_earth-r360 if(dlon_earthmax_rrr) cycle end if ! Extract radial wind data height= thishgt rwnd = thisvr azm_earth = corrected_azimuth if(regional) then if(oneobtest .and. learthrel_rw) then ! for non rotated winds!!! cosazm=cos(azm_earth*deg2rad) sinazm=sin(azm_earth*deg2rad) azm=atan2(sinazm,cosazm)*rad2deg else cosazm_earth=cos(azm_earth*deg2rad) sinazm_earth=sin(azm_earth*deg2rad) call rotate_wind_ll2xy(cosazm_earth,sinazm_earth,cosazm,sinazm,dlon_earth,dlon,dlat) azm=atan2(sinazm,cosazm)*rad2deg end if else azm=azm_earth end if !#################### Data thinning ################### icntpnt=icntpnt+1 ithin=1 !number of obs to keep per grid box if(radar_no_thinning) then ithin=-1 endif if(ithin > 0)then if(zflag == 0)then klon1= int(dlon); klat1= int(dlat) dx = dlon-klon1; dy = dlat-klat1 dx1 = one-dx; dy1 = one-dy w00=dx1*dy1; w10=dx1*dy; w01=dx*dy1; w11=dx*dy klat1=min(max(1,klat1),nlat); klon1=min(max(0,klon1),nlon) if (klon1==0) klon1=nlon klatp1=min(nlat,klat1+1); klonp1=klon1+1 if (klonp1==nlon+1) klonp1=1 do kk=1,nsig hges(kk)=w00*hgtl_full(klat1 ,klon1 ,kk) + & w10*hgtl_full(klatp1,klon1 ,kk) + & w01*hgtl_full(klat1 ,klonp1,kk) + & w11*hgtl_full(klatp1,klonp1,kk) end do sin2 = sin(thislat)*sin(thislat) termg = grav_equator * & ((one+somigliana*sin2)/sqrt(one-eccentricity*eccentricity*sin2)) termr = semi_major_axis /(one + flattening + grav_ratio - & two*flattening*sin2) termrg = (termg/grav)*termr do kk=1,nsig zges(kk) = (termr*hges(kk)) / (termrg-hges(kk)) zl_thin(kk)=zges(kk) end do endif zobs = height ntmp=ndata ! counting moved to map3gridS if (l4dvar) then timedif = zero else timedif=abs(t4dvo-toff) endif crit1 = timedif/r6+half call map3grids(1,zflag,zl_thin,nlevz,dlat_earth,dlon_earth,& zobs,crit1,ndata,iout,icntpnt,iiout,luse, .false., .false.) maxout=max(maxout,iout) maxdata=max(maxdata,ndata) if (.not. luse) then ntdrvr_thin2=ntdrvr_thin2+1 cycle endif if(iiout > 0) isort(iiout)=0 if (ndata > ntmp) then nodata=nodata+1 endif isort(icntpnt)=iout else ndata =ndata+1 nodata=nodata+1 iout=ndata isort(icntpnt)=iout endif !#################### Data thinning ################### if(.not. oneobtest) then iaaa=azm/(r360/(r8*irrr)) iaaa=mod(iaaa,8*irrr) if(iaaa<0) iaaa=iaaa+8*irrr iaaa=iaaa+1 iaaamax=max(iaaamax,iaaa) iaaamin=min(iaaamin,iaaa) end if error = erradar_inflate*thiserr errmax=max(error,errmax) if(thiserr>zero) errmin=min(error,errmin) ! Perform limited qc based on azimuth angle, radial wind ! speed, distance from radar site, elevation of radar, ! height of observation, and observation error good0=.true. if(abs(azm)>r400) then ibadazm=ibadazm+1; good0=.false. end if if(abs(rwnd)>r200) then ibadwnd=ibadwnd+1; good0=.false. end if if(dist>r400) then ibaddist=ibaddist+1; good0=.false. end if if(staheight<-r1000.or.staheight>r50000) then ibadstaheight=ibadstaheight+1; good0=.false. end if if(height<-r1000.or.height>r50000) then ibadheight=ibadheight+1; good0=.false. end if if(heightr6 .or. thiserr<=zero) then ibaderror=ibaderror+1; good0=.false. end if good=.true. if(.not.good0) then notgood0=notgood0+1 cycle end if ! If data is good, load into output array if(good) then usage = zero if(icuse(ikx) < 0)usage=r100 nsuper2_kept=nsuper2_kept+1 cdata(1) = error ! wind obs error (m/s) cdata(2) = dlon ! grid relative longitude cdata(3) = dlat ! grid relative latitude cdata(4) = height ! obs absolute height (m) cdata(5) = rwnd ! wind obs (m/s) cdata(6) = azm*deg2rad ! azimuth angle (radians) cdata(7) = t4dvo+time_offset ! obs time (hour) cdata(8) = ikx ! type cdata(9) = tiltangle ! tilt angle (radians) cdata(10)= staheight ! station elevation (m) cdata(11)= rstation_id ! station id cdata(12)= usage ! usage parameter cdata(13)= idomsfc ! dominate surface type cdata(14)= skint ! skin temperature cdata(15)= ff10 ! 10 meter wind factor cdata(16)= sfcr ! surface roughness cdata(17)=dlon_earth*rad2deg ! earth relative longitude (degrees) cdata(18)=dlat_earth*rad2deg ! earth relative latitude (degrees) cdata(19)=dist ! range from radar in km (used to estimatebeam spread) cdata(20)=zsges ! model elevation at radar site cdata(21)=thiserr cdata(22)=two do j=1,maxdat cdata_all(j,iout)=cdata(j) end do else notgood = notgood + 1 end if end do end do end do close(lnbufr) ! A simple unformatted fortran file should not be mixed with bufr I/O write(6,*)'READ_RADAR_L2RW_NOVADQC: ',trim(outmessage),' reached eof on 2 superob radar file' write(6,*)'READ_RADAR_L2RW_NOVADQC: nsuper2_in,nsuper2_kept=',nsuper2_in,nsuper2_kept write(6,*)'READ_RADAR_L2RW_NOVADQC: # bad winds =',ibadwnd,nobs_badvr,nobs_badsr write(6,*)'READ_RADAR_L2RW_NOVADQC: # num thinned =',kthin,ntdrvr_thin2 write(6,*)'READ_RADAR_L2RW_NOVADQC: timemin,max =',timemin,timemax write(6,*)'READ_RADAR_L2RW_NOVADQC: errmin,max =',errmin,errmax write(6,*)'READ_RADAR_L2RW_NOVADQC: dlatmin,max,dlonmin,max=',dlatmin,dlatmax,dlonmin,dlonmax ! Write observation to scratch file call count_obs(ndata,maxdat,ilat,ilon,cdata_all,nobs) write(lunout) obstype,sis,nreal,nchanl,ilat,ilon write(6,*) shape(cdata_all) write(lunout) ((cdata_all(k,i),k=1,maxdat),i=1,ndata) deallocate(cdata_all) if (radar_sites) deallocate(rtable,rsite,ruse) deallocate(zl_thin) deallocate(isort) return end subroutine read_radar_l2rw !!!!!!!!!!!!!!! End added for l2rw thinning !!!!!!!!!!!!!!! subroutine get_azimuth_corrected(this_stalon,this_stalat, & thislon, thislat, & thisazm, corrected_azm) !=======================================================================! ! Author: G. Zhao ! Follow the alogrithm in subroutine radar_bufr_read_all (in ! read_l2bufr_mod.f90) to get corrected azimuth ! ! Note: all input and output are in unit of degree, not radian. !=======================================================================! use kinds, only: r_kind,r_single,i_kind use constants, only: one,half,two,deg2rad,rad2deg,zero_single implicit none real(r_kind) , intent(in ) :: this_stalon,this_stalat,thislon,thislat,thisazm real(r_kind) , intent(out ) :: corrected_azm ! local variables real(r_kind) rlon0,clat0,slat0,rlonglob,rlatglob,clat1,caz0,saz0,cdlon,sdlon,caz1,saz1 real(r_kind) this_stalatr,thisazimuthr rlon0=this_stalon*deg2rad this_stalatr=this_stalat*deg2rad rlonglob=thislon*deg2rad rlatglob=thislat*deg2rad thisazimuthr=thisazm*deg2rad ! Get corrected azimuth clat0=cos(this_stalatr) slat0=sin(this_stalatr) clat1=cos(rlatglob) caz0=cos(thisazimuthr) saz0=sin(thisazimuthr) cdlon=cos(rlonglob-rlon0) sdlon=sin(rlonglob-rlon0) caz1=clat0*caz0/clat1 saz1=saz0*cdlon-caz0*sdlon*slat0 corrected_azm=atan2(saz1,caz1)*rad2deg return end subroutine get_azimuth_corrected ! ! ! subroutine get_rdr_obshgttilt(thistilt,thisrange,this_stahgt, & thishgt,corrected_tilt) !-------------------------------------------------------------------! ! following the algorithm used in read_l2bufr_mod.f90 to calculate the ! observation height (radar beam height). ! ! 2017-10-19 CAPS(G. Zhao) ! Input: ! thistilt scan tilt of radar beam (elevation angle at radar station, ! degree) ! thisrange radar beam range at observation point (meter) ! this_stahgt height of radar station (meter) ! ! Output: ! thishgt radar beam height at observation point (meter) ! corrected_tilt elevation angle of radar beam at observation point ! (degree, sometimes called corrected tilt angle) ! Due to curvature of radar beam and surface of earth !-------------------------------------------------------------------! ! use kinds, only: r_kind,r_single,i_kind use constants, only: zero,half,one,two,rearth,deg2rad,rad2deg,zero_quad,one_quad implicit none real(r_kind) ,intent(in ) :: thistilt, thisrange, this_stahgt real(r_kind) ,intent( out) :: thishgt real(r_kind) ,intent( out) :: corrected_tilt ! local parameters real(r_kind),parameter:: four_thirds = 4.0_r_kind / 3.0_r_kind real(r_kind),parameter:: r8 = 8.0_r_kind real(r_kind),parameter:: r89_5 = 89.5_r_kind real(r_kind),parameter:: r90 = 90.0_r_kind real(r_kind),parameter:: missval = -999.0_r_kind ! local variables real(r_kind) a43,aactual,b,c,selev0,celev0,epsh,erad,h,ha real(r_kind) rad_per_meter real(r_kind) thistiltr real(r_kind) celev,selev rad_per_meter= one/rearth erad = rearth ! use 4/3rds rule to get elevation of radar beam ! (if local temperature available, then vertical position can be ! estimated with greater accuracy) aactual=erad+this_stahgt a43=four_thirds*aactual thistiltr=thistilt*deg2rad selev0=sin(thistiltr) celev0=cos(thistiltr) b=thisrange*(thisrange+two*aactual*selev0) c=sqrt(aactual*aactual+b) ha=b/(aactual+c) epsh=(thisrange*thisrange-ha*ha)/(r8*aactual) h=ha-epsh thishgt=this_stahgt+h if(thishgt < zero)then ! don't use observation if it is likely to be below sealevel thishgt=missval end if ! Get corrected tilt angle (i.e., elevation angle at obs location) celev=celev0 selev=selev0 if(thisrange>=one) then celev=a43*celev0/(a43+h) selev=(thisrange*thisrange+h*h+two*a43*h)/(two*thisrange*(a43+h)) end if corrected_tilt=atan2(selev,celev)*rad2deg return end subroutine get_rdr_obshgttilt !===============================================================================! ! subroutines from CAPS ARPS package ! adas/radarlib3d.f90 ! including: ! subroutine beamelv dhdrange !===============================================================================! ! subroutine beamelv(height,sfcrng,elvang,range) ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Calculate the elevation angle (elvang) and the along ! ray-path distance (range) of a radar beam ! crossing through the given height and along-ground ! distance. ! ! This method assumes dn/dh is constant such that the ! beam curves with a radius of 4/3 of the earth's radius. ! This is dervied from Eq. 2.28 of Doviak and Zrnic', ! Doppler Radar and Weather Observations, 1st Ed. ! !----------------------------------------------------------------------- ! ! AUTHOR: Keith Brewster ! 10/10/95 ! ! MODIFICATION HISTORY: ! !----------------------------------------------------------------------- ! ! INPUT: ! height Height (meters) of beam above ground. ! sfcrng Distance (meters) of point along ground from radar. ! ! OUTPUT ! elvang Elevation angle (degrees) of radar beam ! range Distance (meters) along radar beam from radar ! !----------------------------------------------------------------------- ! ! !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- ! use kinds, only: r_kind,r_single,i_kind,r_double implicit none real(r_kind), intent(in ) :: height real(r_kind), intent(in ) :: sfcrng real(r_kind), intent( out) :: elvang real(r_kind), intent( out) :: range ! real(r_double) :: eradius,frthrde,rad2deg parameter (eradius=6371000._r_double, & frthrde=(4._r_double*eradius/3._r_double), & rad2deg=(180._r_double/3.14592654_r_double)) ! real(r_double):: elvrad,hgtdb,rngdb,drange ! if (sfcrng > 0._r_kind) then hgtdb=frthrde+dble(height) rngdb=dble(sfcrng)/frthrde elvrad = atan((hgtdb*cos(rngdb) - frthrde)/(hgtdb * sin(rngdb))) drange = (hgtdb*sin(rngdb))/cos(elvrad) elvang=rad2deg*elvrad range=drange else elvang=90._r_kind range=height end if return end subroutine beamelv ! subroutine dhdrange(elvang,range,dhdr) ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Calculate the local change in height of the radar ! beam with respect to a change in range. Due to ! curvature of the beam and the earth's surface this is ! generally different what would be calculated from the ! elevation angle measured at the radar. This derivative ! is needed for finding 3-d velocities from radial winds ! and accounting for terminal velocity of precipitation. ! ! This formulation, consistent with subroutine beamhgt, ! assumes a 4/3 earth radius beam curvature. This formula ! is obtained by differentiating Eq 2.28 of Doviak and ! Zrnic', Doppler Radar and Weather Observations, 1st Ed. ! !----------------------------------------------------------------------- ! ! AUTHOR: Keith Brewster ! 06/22/95 ! ! MODIFICATION HISTORY: ! ! !----------------------------------------------------------------------- ! ! INPUT: ! ! elvang Elevation angle (degrees) of radar beam ! range Distance (meters) along radar beam from radar ! ! OUTPUT: ! dhdr Change in height per change in range (non-dimensional) ! ! !----------------------------------------------------------------------- ! ! !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- ! use kinds, only: r_kind,r_single,i_kind,r_double implicit none real(r_kind), intent(in ) :: range real(r_kind), intent(in ) :: elvang real(r_kind), intent( out) :: dhdr ! real(r_double) :: eradius,frthrde,eighthre,fthsq,deg2rad parameter (eradius=6371000._r_double, & frthrde=(4._r_double*eradius/3._r_double), & eighthre=(8._r_double*eradius/3._r_double), & fthsq=(frthrde*frthrde), & deg2rad=(3.14592654_r_double/180._r_double)) ! real(r_double) :: sinelv,dhdrdb,drange ! drange=dble(range) sinelv=sin(deg2rad*dble(elvang)) dhdrdb = (drange+frthrde*sinelv)/ & sqrt(drange*drange + fthsq + eighthre*drange*sinelv) dhdr = dhdrdb ! return end subroutine dhdrange