subroutine read_radar_wind_ascii(nread,ndata,nodata,infile,lunout,obstype,sis,hgtl_full,nobs) !$$$ subprogram documentation block ! . . . . ! subprogram: read_dbz read level2 raw QC'd radial velocity data ! ! prgmmr: carley org: np22 date: 2011-05-24 ! ! abstract: Reads and processes level 2 horizontal radial velocity (m/s) by ! radar site. Data are on radar scan surafces. Processing includes ! finding the lat/lon and height of each observation. . ! ! program history log: ! 2011-08-12 carley - fix ob error to 2 m/s ! 2011-08-23 carley - use deter_sfc_mod ! 2011-12-08 carley - add wind rotation (earth to grid) ! ! input argument list: ! infile - file from which to read data ! lunout - unit to which to write data for further processing ! obstype - observation type to process ! ! output argument list: ! nread - number of radar reflectivity observations read ! ndata - number of radar reflectivity observations retained for further processing ! nodata - number of radar reflectivity observations retained for further processing ! sis - satellite/instrument/sensor indicator ! ! Variable Definitions: ! ! a43 - real - (4/3)*(earth radius) ! a,b,c,ha,epsh,h,aactual - real - used in computing radar observation height ! cdata_all - real - dim(maxdat,maxobs) - array holding all data for assimilation ! celev0,selev0 - real- cos and sin of elevation angle (raw) ! celev,selev - real - corrected cos and sin of elevation angle ! clat0 - real - cos of radar station latitude ! cstaid - char - radar station ide ! dbzerr - real - observation error (obtained from convinfo - dBZ) ! dlat - real - grid relative latitude of observation (grid units) ! dlon - real - grid relative longitude of observation (grid units) ! gamma - real - used in finding observation latlon ! lunrad - int - unit number for reading radar data from file ! maxobs - int - max number of obs converted to no precip observations ! num_m2nopcp -int - number of missing obs ! num_missing - int - number of missing observations ! num_noise - int - number of rejected noise observations ! num_nopcp - int - number of noise obs converted to no precip observations ! numbadtime - int - number of elevations outside time window ! num_badtilt - int - number of elevations outside specified interval ! num_badrange - int - number of obs outside specified range distance ! obdate - int - dim(5) - yyyy,mm,dd,hh,minmin of observation ! outside - logical - if observations are outside the domain -> true ! radartwindow - real - time window for radar observations (minutes) ! rlatglob - real - earth relative latitude of observation (radians) ! rlatloc - real - latitude of observation on radar-relative projection ! rlonglob - real - earth relative longitude of observation (radians) ! rlonloc - real - longitude of observation on radar-relative projection ! rlon0 - real - radar station longitude (radians) ! rmins_an - real - analysis time from reference date (minutes) ! rmins_ob - real - observation time from reference date (minutes) ! rstation_id - real - radar station id ! slat0 - real - sin of radar station latitude ! thisazimuthr - real - 90deg minues the actual azimuth and converted to radians ! thiserr - real - observation error ! thislat - real - latitude of observation ! thislon - real - longitude of observation ! thisrange - real - range of observation from radar ! thishgt - real - observation height ! this_stahgt - real - radar station height (meters about sea level) ! this_staid - char - radar station id ! thistilt - real - radar tilt angle (degrees) ! thistiltr - real- radar tilt angle (radians) ! timeb - real - obs time (analyis relative minutes) ! ! ! ! Derived data types ! ! radar - derived data type for containing volume scan information ! nelv- int - number of elevation angles ! radid - char*4 - radar ID (e.g. KAMA) ! vcpnum - int - volume coverage pattern number ! year - int - UTC ! day - int - UTC ! month - int - UTC ! hour - in - UTC ! minute - int - UTC ! second - int - UTC ! radhgt - real - elevation of the radar above sea level in meters (I believe ! this includes the height of the antenna as well) ! radlat - real - latitude location of the radar ! radlon - real - longitude location of the radar ! fstgatdis - real - first gate distance (meters) ! gatewidth - real - gate width (meters) ! elev_angle - real - radar elevation angle (degrees) ! num_beam - int - number of beams ! num_gate - int - number of gates ! nyq_vel - real - nyquist velocity ! azim - real - azimuth angles ! field - real - radar data variable (reflectivity or velocity) ! ! Defined radar types: ! strct_in_vel - radar - contains volume scan information related to ! radial velocity ! strct_in_dbz - radar - contains volume scan information related to ! radar reflectivity ! strct_in_rawvel - radar - contains volume scan information related to ! raw radial velocity ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block use kinds, only: r_kind,r_double,i_kind use constants, only: zero,half,one,two,deg2rad,rearth,rad2deg, & one_tenth,r1000,r60,r60inv,r100,r400,grav_equator, & eccentricity,somigliana,grav_ratio,grav,semi_major_axis,flattening use gridmod, only: regional,tll2xy,rotate_wind_ll2xy,nsig,nlat,nlon use obsmod, only: iadate, & mintiltvr,maxtiltvr,minobrangevr,maxobrangevr, rmesh_vr,zmesh_vr,& doradaroneob,oneoblat,oneoblon,oneobheight,oneobradid use obsmod,only: radar_no_thinning use gsi_4dvar, only: l4dvar,time_4dvar use convinfo, only: nconvtype,ctwind,icuse,ioctype use convthin, only: make3grids,map3grids,del3grids,use_all use read_l2bufr_mod, only: invtllv use qcmod, only: erradar_inflate use deter_sfc_mod, only: deter_sfc2,deter_zsfc_model use mpimod, only: npe implicit none ! Declare passed variables character(len=*),intent(in ) :: obstype,infile character(len=*),intent(in ) :: sis integer(i_kind) ,intent(in ) :: lunout integer(i_kind) ,intent(inout) :: nread,ndata,nodata real(r_kind),dimension(nlat,nlon,nsig),intent(in):: hgtl_full integer(i_kind),dimension(npe) ,intent(inout) :: nobs ! Declare 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:: r6 = 6.0_r_kind real(r_kind),parameter:: r360=360.0_r_kind integer(i_kind),parameter:: maxdat=22 ! Used in generating cdata array !--Derived data type declaration type :: radar character(4) :: radid integer(i_kind) :: vcpnum integer(i_kind) :: year integer(i_kind) :: month integer(i_kind) :: day integer(i_kind) :: hour integer(i_kind) :: minute integer(i_kind) :: second real(r_kind) :: radlat real(r_kind) :: radlon real(r_kind) :: radhgt real(r_kind) :: fstgatdis real(r_kind) :: gateWidth real(r_kind) :: elev_angle integer(i_kind) :: num_beam integer(i_kind) :: num_gate real(r_kind) :: nyq_vel real(r_kind),allocatable :: azim(:) !has dimension (num_beam) real(r_kind),allocatable :: field(:,:) !has dimension (num_gate,num_beam) end type radar !--Counters for diagnostics integer(i_kind) :: num_missing=0,numbadtime=0, & !counts num_badtilt=0,num_badrange=0, & ibadazm=0 integer(i_kind) :: ithin,zflag,nlevz,icntpnt,klon1,klat1,kk,klatp1,klonp1 real(r_kind) :: rmesh,xmesh,zmesh,dx,dy,dx1,dy1,w00,w01,w10,w11 real(r_kind), allocatable, dimension(:) :: zl_thin real(r_kind),dimension(nsig):: hges,zges real(r_kind) sin2,termg,termr,termrg,zobs,height integer(i_kind) ntmp,iout,iiout,ntdrvr_thin2 real(r_kind) crit1,timedif real(r_kind),parameter:: r16000 = 16000.0_r_kind logical :: luse integer(i_kind) maxout,maxdata integer(i_kind),allocatable,dimension(:):: isort !--General declarations integer(i_kind) :: ierror,lunrad,i,j,k,v,na,nb,nelv,nvol, & ikx,mins_an,mins_ob integer(i_kind) :: maxobs,nchanl,ilat,ilon,idomsfc integer(i_kind),dimension(5) :: obdate real(r_kind) :: b,c,ha,epsh,h,aactual,a43,thistilt,ff10,sfcr,skint,zsges, & radar_lon,radar_lat,dlon_radar,dlat_radar,errmax,errmin,error real(r_kind) :: thistiltr,selev0,celev0,thisrange,this_stahgt,thishgt real(r_kind) :: celev,selev,gamma,thisazimuthr,rlon0,t4dv, & clat0,slat0,dlat,dlon,thiserr,thislon,thislat, & rlonloc,rlatloc,rlonglob,rlatglob,timeb,rad_per_meter real(r_kind) :: azm,cosazm_earth,sinazm_earth,cosazm,sinazm real(r_kind) :: radartwindow real(r_kind) :: rmins_an,rmins_ob 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) logical :: outside type(radar),allocatable :: strct_in_vel(:,:) real(r_kind) :: mintilt,maxtilt,maxobrange,minobrange integer(i_kind) :: thin_freq=1 mintilt=mintiltvr maxtilt=maxtiltvr minobrange=minobrangevr maxobrange=maxobrangevr !-Check if radial velocity is in the convinfo file and extract necessary attributes ithin=1 !number of obs to keep per grid box if(radar_no_thinning) then ithin=-1 endif errmax=-huge(errmax) errmin=huge(errmin) ikx=0 do i=1,nconvtype if(trim(obstype) == trim(ioctype(i)) .and. abs(icuse(i))== 1) then ikx=i radartwindow=ctwind(ikx)*r60 !Time window units converted to minutes ! (default setting for dbz within convinfo is 0.05 hours) thiserr= 2_r_kind !1.75_r_kind !2_r_kind !Ob error (m/s) to use for radial velocity exit !Exit loop when finished with initial convinfo fields else if ( i==nconvtype ) then write(6,*) 'READ_RADAR_WIND_ASCII: ERROR - OBSERVATION TYPE IS NOT PRESENT IN CONVINFO OR USE FLAG IS ZERO' write(6,*) 'READ_RADAR_WIND_ASCII: ABORTING read_radar_wind_ascii.f90 - NO VELOCITY OBS READ!' return endif end do if (minobrange >= maxobrange) then write(6,*) 'MININMUM OB RANGE >= MAXIMUM OB RANGE FOR READING RADIAL VELOCITY - PROGRAM STOPPING FROM READ_RADAR_WIND_ASCII.F90' call stop2(400) end if !-next three values are dummy values for now nchanl=0 ilon=2 ilat=3 maxobs=50000000 !value taken from read_radar.f90 !--Allocate cdata_all array allocate(cdata_all(maxdat,maxobs),isort(maxobs)) rmesh=rmesh_vr zmesh=zmesh_vr maxout=0 maxdata=0 isort=0 ntdrvr_thin2=0 icntpnt=0 zflag=0 use_all=.true. if (ithin > 0) then write(6,*)'READ_RADAR: ithin,rmesh :',ithin,rmesh use_all=.false. if(zflag == 0)then nlevz=nsig else nlevz=r16000/zmesh endif xmesh=rmesh call make3grids(xmesh,nlevz) allocate(zl_thin(nlevz)) if (zflag == 1) then do k=1,nlevz zl_thin(k)=k*zmesh enddo endif write(6,*)'READ_RADAR: xmesh, zflag, nlevz =', xmesh, zflag, nlevz endif lunrad=31 open(lunrad,file=trim(infile),status='old',action='read', & iostat=ierror,form='formatted') fileopen: if (ierror == 0) then read(lunrad,'(2i8)') nelv,nvol !read number of elevations and number of volumes !*************************IMPORTANT***************************! ! ! ! All data = 999.0 correspond to missing or bad data ! ! ! !*************************************************************! !------Begin processing--------------------------! !-Obtain analysis time in minutes since reference date call w3fs21(iadate,mins_an) !mins_an -integer number of mins snce 01/01/1978 rmins_an=mins_an !convert to real number volumes: do v=1,nvol read(lunrad,'(i8)') nelv allocate(strct_in_vel(1,nelv)) tilts: do k=1,nelv read(lunrad,'(a4)') strct_in_vel(1,k)%radid read(lunrad,'(i8)') strct_in_vel(1,k)%vcpnum read(lunrad,'(6i8)') strct_in_vel(1,k)%year & ,strct_in_vel(1,k)%month & ,strct_in_vel(1,k)%day & ,strct_in_vel(1,k)%hour & ,strct_in_vel(1,k)%minute & ,strct_in_vel(1,k)%second read(lunrad,'(2f10.3,f10.1)') strct_in_vel(1,k)%radlat & ,strct_in_vel(1,k)%radlon & ,strct_in_vel(1,k)%radhgt read(lunrad,'(2f8.1)') strct_in_vel(1,k)%fstgatdis & ,strct_in_vel(1,k)%gateWidth read(lunrad,'(f8.3)') strct_in_vel(1,k)%elev_angle read(lunrad,'(2i8)') strct_in_vel(1,k)%num_beam & ,strct_in_vel(1,k)%num_gate na=strct_in_vel(1,k)%num_beam nb=strct_in_vel(1,k)%num_gate !******allocate arrays within radar data type**********! allocate(strct_in_vel(1,k)%azim(na)) allocate(strct_in_vel(1,k)%field(nb,na)) !******************************************************! read(lunrad,'(f8.3)') strct_in_vel(1,k)%nyq_vel read(lunrad,'(15f6.1)') (strct_in_vel(1,k)%azim(j),j=1,na) read(lunrad,'(20f6.1)') ((strct_in_vel(1,k)%field(i,j),i=1,nb),j=1,na) obdate(1)=strct_in_vel(1,k)%year obdate(2)=strct_in_vel(1,k)%month obdate(3)=strct_in_vel(1,k)%day obdate(4)=strct_in_vel(1,k)%hour obdate(5)=strct_in_vel(1,k)%minute call w3fs21(obdate,mins_ob) !mins_ob -integer number of mins snce 01/01/1978 rmins_ob=mins_ob !convert to real number rmins_ob=rmins_ob+(strct_in_vel(1,k)%second*r60inv) !convert seconds to minutes and add to ob time !-Comparison is done in units of minutes timeb = rmins_ob-rmins_an if(doradaroneob .and. (oneobradid /= strct_in_vel(1,k)%radid)) cycle tilts if(abs(timeb) > abs(radartwindow)) then numbadtime=numbadtime+1 cycle tilts !If not in time window, cycle the loop end if !--Time window check complete--! thistilt=strct_in_vel(1,k)%elev_angle if (thistilt <= maxtilt .and. thistilt >= mintilt) then gates: do i=1,strct_in_vel(1,k)%num_gate,thin_freq thisrange=strct_in_vel(1,k)%fstgatdis + float(i-1)*strct_in_vel(1,k)%gateWidth !-Check to make sure observations are within specified range if (thisrange <= maxobrange .and. thisrange >= minobrange) then azms: do j=1,strct_in_vel(1,k)%num_beam !-Check to see if this is a missing observation) nread=nread+1 if ( strct_in_vel(1,k)%field(i,j) >= 999.0_r_kind ) then num_missing=num_missing+1 cycle azms !No reason to process the ob if it is missing end if !--Find observation height using method from read_l2bufr_mod.f90 this_stahgt=strct_in_vel(1,k)%radhgt aactual=rearth+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 height=thishgt !--Find observation location using method from read_l2bufr_mod.f90 !-Get corrected tilt angle celev=celev0 selev=selev0 celev=a43*celev0/(a43+h) selev=(thisrange*thisrange+h*h+two*a43*h)/(two*thisrange*(a43+h)) gamma=half*thisrange*(celev0+celev) !-Get earth lat lon of observation rlon0=deg2rad*strct_in_vel(1,k)%radlon clat0=cos(deg2rad*strct_in_vel(1,k)%radlat) slat0=sin(deg2rad*strct_in_vel(1,k)%radlat) thisazimuthr=(90.0_r_kind-strct_in_vel(1,k)%azim(j))*deg2rad !Storing as 90-azm to ! be consistent with ! read_l2bufr_mod.f90 rad_per_meter=one/rearth 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 if(doradaroneob) then thislat=oneoblat thislon=oneoblon thishgt=oneobheight endif if(thislon>=r360) thislon=thislon-r360 if(thislonzero) errmin=min(error,errmin) if(abs(azm)>r400) then ibadazm=ibadazm+1 cycle azms end if this_staid=strct_in_vel(1,k)%radid !Via equivalence in declaration, value is propagated ! to rstation_id used below. ! Get model terrain at radar station location ! If radar station is outside of grid, does not mean the ! radar obs are outside the grid - therefore no need to ! cycle azms. radar_lon=deg2rad*strct_in_vel(1,k)%radlon radar_lat=deg2rad*strct_in_vel(1,k)%radlat call tll2xy(radar_lon,radar_lat,dlon_radar,dlat_radar,outside) call deter_zsfc_model(dlat_radar,dlon_radar,zsges) ! Determines land surface type based on surrounding land ! surface types t4dv=timeb*r60inv call deter_sfc2(thislat,thislon,t4dv,idomsfc,skint,ff10,sfcr) !#################### 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(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(t4dv-toff) timedif=abs(t4dv) !don't know about this endif crit1 = timedif/r6+half call map3grids(1,zflag,zl_thin,nlevz,thislat,thislon,& 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 cdata_all(1,iout) = error ! wind obs error (m/s) cdata_all(2,iout) = dlon ! grid relative longitude cdata_all(3,iout) = dlat ! grid relative latitude cdata_all(4,iout) = thishgt ! obs absolute height (m) cdata_all(5,iout) = strct_in_vel(1,k)%field(i,j) ! wind obs (m/s) cdata_all(6,iout) = azm ! azimuth angle (radians) cdata_all(7,iout) = t4dv ! obs time (hour) - analysis relative cdata_all(8,iout) = ikx ! type cdata_all(9,iout) = thistiltr ! tilt angle (radians) cdata_all(10,iout)= this_stahgt ! station elevation (m) cdata_all(11,iout)= rstation_id ! station id cdata_all(12,iout)= icuse(ikx) ! usage parameter cdata_all(13,iout)= idomsfc ! dominate surface type cdata_all(14,iout)= skint ! skin temperature cdata_all(15,iout)= ff10 ! 10 meter wind factor cdata_all(16,iout)= sfcr ! surface roughness cdata_all(17,iout)=thislon*rad2deg ! earth relative longitude (degrees) cdata_all(18,iout)=thislat*rad2deg ! earth relative latitude (degrees) cdata_all(19,iout)=thisrange/1000_r_kind ! range from radar in km (used to estimate beam spread) cdata_all(20,iout)=zsges ! model elevation at radar site cdata_all(21,iout)=thiserr cdata_all(22,iout)=two ! Level 2 data if(doradaroneob .and. (cdata_all(5,iout) > -99_r_kind) ) exit volumes end do azms !j else num_badrange=num_badrange+1 !If outside acceptable range, increment end if !Range check end do gates !i else num_badtilt=num_badtilt+1 !If outside acceptable tilts, increment end if !Tilt check end do tilts !k do k=1,nelv deallocate(strct_in_vel(1,k)%azim) deallocate(strct_in_vel(1,k)%field) enddo deallocate(strct_in_vel) end do volumes !v close(lunrad) !modified to do one scan at a time if (.not. use_all) then deallocate(zl_thin) call del3grids endif !end modified for thinning !---all looping done now print diagnostic output write(6,*)'READ_RADAR_WIND_ASCII: Reached eof on radar wind ascii file' write(6,*)'READ_RADAR_WIND_ASCII: # volumes in input file =',nvol write(6,*)'READ_RADAR_WIND_ASCII: # elevations per volume =',nelv write(6,*)'READ_RADAR_WIND_ASCII: # elevations outside time window =',numbadtime write(6,*)'READ_RADAR_WIND_ASCII: # of missing data =',num_missing write(6,*)'READ_RADAR_WIND_ASCII: # outside specif. range =',num_badrange write(6,*)'READ_RADAR_WIND_ASCII: # outside specif. tilts =',num_badtilt write(6,*)'READ_RADAR_WIND_ASCII: # bad azimuths =',ibadazm !---Write observation to scratch file---! call count_obs(ndata,maxdat,ilat,ilon,cdata_all,nobs) write(lunout) obstype,sis,maxdat,nchanl,ilat,ilon write(lunout) ((cdata_all(k,i),k=1,maxdat),i=1,ndata) !---------------DEALLOCATE ARRAYS-------------! deallocate(cdata_all) else !fileopen write(6,*) 'READ_RADAR_WIND_ASCII: ERROR OPENING RADIAL VELOCITY FILE: ',trim(infile),' IOSTAT ERROR: ',ierror, ' SKIPPING...' end if fileopen end subroutine read_radar_wind_ascii