subroutine read_dbz_mrms_netcdf(nread,ndata,nodata,infile,obstype,lunout,sis,nobs) !$$$ subprogram documentation block ! . . . . ! subprogram: read_dbz read level2 raw QC'd radar reflectivity files ! ! prgmmr: carley org: np22 date: 2011-04-04 ! ! abstract: Reads and processes level 2 horizontal radar reflectivity (dBZ) by ! radar site. Data are on radar scan surafces. Also reads, but does ! not process unfolded radial velocities. Processing includes ! finding the lat/lon and height of each observation. ! This formulation is not outfitted for 4dvar, but will ! work with 3dvar and hybrid ensemble. ! ! program history log: ! 2011-08-12 carley - Fix dBZ oberror to be 3dBZ and add optional ! upper bound limit to observed dBZ to account ! for representativeness error. ! 2011-12-08 carley - Fix dBZ oberror to 5 dBZ ! ! 2015 Lei -- modify from read_dbz to read_dbz_mrms_netcdf ! 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 netcdf use kinds, only: r_kind,r_double,i_kind,r_single use constants, only: zero,half,one,two,deg2rad,rearth,rad2deg, & one_tenth,r1000,r60,r60inv,r100,r400 use gridmod, only: tll2xy use obsmod, only: iadate use convinfo, only: nconvtype,ctwind,icuse,ioctype use mpimod, only: npe use read_l2bufr_mod, only : invtllv 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 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:: r360=360.0_r_kind integer(i_kind),parameter:: maxdat=17 ! Used in generating cdata array integer (i_kind):: iyear,imon,iday,ihour,imin,isec !--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,num_nopcp=0, & !counts numbadtime=0,num_badtilt=0, & num_badrange=0,num_m2nopcp=0, & num_noise=0,num_limmax=0 ,num_limmin=0 !--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,scount integer(i_kind),dimension(5) :: obdate real(r_kind) :: b,c,ha,epsh,h,aactual,a43,thistilt real(r_kind) :: thistiltr,selev0,celev0,thisrange,this_stahgt,thishgt real(r_kind) :: celev,selev,gamma,thisazimuthr,rlon0, & clat0,slat0,dlat,dlon,thiserr,thislon,thislat, & rlonloc,rlatloc,rlonglob,rlatglob,timeb,rad_per_meter real(r_kind) :: radartwindow real(r_kind) :: dbzerr,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_dbz(:,:) !---------SETTINGS FOR FUTURE NAMELIST---------! integer(i_kind) :: maxobrange=999000 ! Range (m) *within* which to use observations - obs *outside* this range are not used integer(i_kind) :: minobrange=-999 ! Range (m) *outside* of which to use observatons - obs *inside* this range are not used real(r_kind) :: mintilt=0.0_r_kind ! Only use tilt(elevation) angles (deg) >= this number real(r_kind) :: maxtilt=20.0_r_kind ! Do no use tilt(elevation) angles (deg) >= this number logical :: missing_to_nopcp=.false. ! Set missing observations to 'no precipitation' observations -> dbznoise (See Aksoy et al. 2009, MWR) real(r_kind) :: dbznoise=2.0_r_kind ! dBZ obs must be >= dbznoise for assimilation logical :: l_limmax=.true. ! If true, observations > 60 dBZ are limited to be 60 dBZ. This is logical :: l_limmin=.true. ! If true, observations <0 dBZ are limited to be 0 dBZ. This is character (len=4) :: radarsite_nc character (len=256) vcpstr_nc !following the treatment on the precision issue for netcdf like in !wrf_netcdf_interface.F90 integer(i_kind) :: ncid,ierr,dimid1,dimid2 integer(i_kind) :: varid1,varid2,varid3,varid4,varid6 integer(i_kind) :: numazim_nc,numgate_nc,vcp_nc real(r_single) :: elev_nc,firstgate_nc,lat_nc,lon_nc,height_nc real(r_single), allocatable :: azimuth_nc(:),beamwidth_nc(:),azimspacing_nc(:),gatewidth_nc(:) real(r_single), allocatable :: nyquist_nc(:),obdata_nc(:,:) real(r_single) nyquist_default_nc parameter(nyquist_default_nc=50.0_r_kind) !clg ! ! due to representativeness error associated with the model !----------------------------------------------! !--------------------------------------------------------------------------------------! ! END OF ALL DECLARATIONS ! !--------------------------------------------------------------------------------------! !-Check if reflectivity is in the convinfo file and extract necessary attributes scount=0 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) dbzerr=5_r_kind !Ob error (dB) to use for radar reflectivity factor exit !Exit loop when finished with initial convinfo fields else if ( i==nconvtype ) then write(6,*) 'READ_dBZ: ERROR - OBSERVATION TYPE IS NOT PRESENT IN CONVINFO OR USE FLAG IS ZERO' write(6,*) 'READ_dBZ: ABORTTING read_dbz.f90 - NO REFLECTIVITY OBS READ!' return endif end do if (minobrange >= maxobrange) then write(6,*) 'MININMUM OB RANGE >= MAXIMUM OB RANGE FOR READING dBZ - PROGRAM STOPPING FROM READ_DBZ.F90' call stop2(400) end if !-next three values are dummy values for now nchanl=0 ilon=2 ilat=3 maxobs=2000000 !value taken from read_radar.f90 !--Allocate cdata_all array allocate(cdata_all(maxdat,maxobs)) lunrad=31 ! get the time from the file name !read the QC radar data in titls in netcdf format from K. Cooper. nvol=1 nelv=1 v=1;k=1 allocate(strct_in_dbz(nvol,nelv)) !!READ RADAR DATA ierr = NF90_OPEN(trim(infile),0,ncid) if (ierr /= nf90_noerr) call handle_err(ierr,"open") ierr = NF90_INQ_DIMID(ncid,'Azimuth',dimid1) if (ierr /= nf90_noerr) call handle_err(ierr,"Azimuth") ierr = NF90_INQ_DIMID(ncid,'Gate',dimid2) if (ierr /= nf90_noerr) call handle_err(ierr,"Gate") ierr = NF90_INQ_VARID(ncid,'Azimuth',varid1) if (ierr /= nf90_noerr) call handle_err(ierr,"Azimuth") ierr = NF90_INQ_VARID(ncid,'BeamWidth',varid2) if (ierr /= nf90_noerr) call handle_err(ierr,"BeamWidth") ierr = NF90_INQ_VARID(ncid,'AzimuthalSpacing',varid3) if (ierr /= nf90_noerr) call handle_err(ierr,"azimuthalspacing") ierr = NF90_INQ_VARID(ncid,'GateWidth',varid4) if (ierr /= nf90_noerr) call handle_err(ierr,"gatewidth") ierr = NF90_INQ_VARID(ncid,'ReflectivityQC',varid6) if (ierr /= nf90_noerr) call handle_err(ierr,"ReflectivityQC") ierr = nf90_inquire_dimension(ncid, dimid1, len = numazim_nc) if (ierr /= nf90_noerr) call handle_err(ierr,"numazim data") ierr = nf90_inquire_dimension(ncid, dimid2, len = numgate_nc) if (ierr /= nf90_noerr) call handle_err(ierr,"numgate data") ierr = NF90_GET_ATT(ncid,nf90_global,'Elevation',elev_nc) if (ierr /= nf90_noerr) call handle_err(ierr,"get elev") ierr = NF90_GET_ATT(ncid,nf90_global,'RangeToFirstGate',firstgate_nc) if (ierr /= nf90_noerr) call handle_err(ierr,"get firstgate") ierr = NF90_GET_ATT(ncid,nf90_global,'Latitude',lat_nc) if (ierr /= nf90_noerr) call handle_err(ierr,"get lat") ierr = NF90_GET_ATT(ncid,nf90_global,'Longitude',lon_nc) if (ierr /= nf90_noerr) call handle_err(ierr,"get lon") ierr = NF90_GET_ATT(ncid,nf90_global,'radarName-value',radarsite_nc) if (ierr /= nf90_noerr) call handle_err(ierr,"radarsite") ierr = NF90_GET_ATT(ncid,nf90_global,'vcp-value',vcpstr_nc) if (ierr /= nf90_noerr) call handle_err(ierr,"vcp") read(vcpstr_nc,*) vcp_nc ierr = NF90_GET_ATT(ncid,nf90_global,'Height',height_nc) if (ierr /= nf90_noerr) call handle_err(ierr,"height") !reverse order of dimensions as stated in ncdump: allocate(azimuth_nc(numazim_nc),beamwidth_nc(numazim_nc),azimspacing_nc(numazim_nc),gatewidth_nc(numazim_nc)) allocate(nyquist_nc(numazim_nc),obdata_nc(numgate_nc,numazim_nc)) ierr = NF90_GET_VAR(ncid,varid1,azimuth_nc) if (ierr /= nf90_noerr) call handle_err(ierr,"azimuth data") ierr = NF90_GET_VAR(ncid,varid2,beamwidth_nc) if (ierr /= nf90_noerr) call handle_err(ierr,"beamwidth data") ierr = NF90_GET_VAR(ncid,varid3,azimspacing_nc) if (ierr /= nf90_noerr) call handle_err(ierr,"azimspacing data") ierr = NF90_GET_VAR(ncid,varid4,gatewidth_nc) if (ierr /= nf90_noerr) call handle_err(ierr,"gatewidth data") ierr = NF90_GET_VAR(ncid,varid6,obdata_nc) if (ierr /= nf90_noerr) call handle_err(ierr,"obdata data") ierr = NF90_CLOSE(ncid) if (ierr /= nf90_noerr) call handle_err(ierr,"close") do i=1,numazim_nc if ( (beamwidth_nc(i) /= beamwidth_nc(1)) .or. (gatewidth_nc(i) /= gatewidth_nc(1)) )then print *, "stopping: non-uniform scan" endif enddo read(infile(21:24),'(I4.4)')iyear read(infile(25:26),'(I2.2)')imon read(infile(27:28),'(I2.2)')iday read(infile(30:31),'(I2.2)')ihour read(infile(32:33),'(I2.2)')imin read(infile(34:35),'(I2.2)')isec do i=1,numgate_nc do j=1,numazim_nc if(obdata_nc(i,j) <= -999_r_kind) obdata_nc(i,j)=-999_r_kind enddo enddo !transform the read-in ob to the intermidate obs variables( radar obs to be used in GSI strct_in_dbz(v,k)%radid=radarsite_nc strct_in_dbz(v,k)%vcpnum=vcp_nc strct_in_dbz(v,k)%year=iyear ! to be defind from infile name strct_in_dbz(v,k)%month=imon strct_in_dbz(v,k)%day=iday strct_in_dbz(v,k)%hour=ihour strct_in_dbz(v,k)%minute=imin strct_in_dbz(v,k)%second=isec strct_in_dbz(v,k)%radlat=lat_nc strct_in_dbz(v,k)%radlon=lon_nc strct_in_dbz(v,k)%radhgt=height_nc strct_in_dbz(v,k)%fstgatdis =firstgate_nc strct_in_dbz(v,k)%gateWidth=gatewidth_nc(1) ! always the same ??) strct_in_dbz(v,k)%elev_angle=elev_nc strct_in_dbz(v,k)%num_beam=numazim_nc strct_in_dbz(v,k)%num_gate=numgate_nc na=strct_in_dbz(v,k)%num_beam nb=strct_in_dbz(v,k)%num_gate !******allocate arrays within radar data type**********! allocate(strct_in_dbz(v,k)%azim(na)) allocate(strct_in_dbz(v,k)%field(nb,na)) !******************************************************! strct_in_dbz(v,k)%azim(:)=azimuth_nc(:) strct_in_dbz(v,k)%field(:,:)=obdata_nc(:,:) ierror=0 fileopen: if (ierror == 0) then !Check to make sure file is open - will also fail if file does not exist. Closing endif at end of subroutine. !*************************IMPORTANT***************************! ! ! ! All data = 999.0 correspond to missing or bad data ! ! ! !*************************************************************! !------Begin processing--------------------------! !-Obtain analysis time in minutes since reference date if(ndata/=0) then ! for further thinking write(6,*)'ndata is not 0 in read_dbz_netcdf, its impact needs to be considered ,stop' endif call w3fs21(iadate,mins_an) !mins_an -integer number of mins snce 01/01/1978 ! w3movedat to help get a date from the time difference rmins_an=mins_an !convert to real number volumes: do v=1,nvol tilts: do k=1,nelv !--Check if observation fits within specified time window--! !-Find reference time of observation obdate(1)=strct_in_dbz(v,k)%year obdate(2)=strct_in_dbz(v,k)%month obdate(3)=strct_in_dbz(v,k)%day obdate(4)=strct_in_dbz(v,k)%hour obdate(5)=strct_in_dbz(v,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_dbz(v,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(abs(timeb) > 100_r_kind) cycle write(6,*) 'Processing obdate:',obdate,strct_in_dbz(v,k)%second !--Time window check complete--! thistilt=strct_in_dbz(v,k)%elev_angle if (thistilt <= maxtilt .and. thistilt >= mintilt) then gates: do i=1,strct_in_dbz(v,k)%num_gate thisrange=strct_in_dbz(v,k)%fstgatdis + float(i-1)*strct_in_dbz(v,k)%gateWidth !-Check to make sure observations are within specified range if (thisrange <= maxobrange .and. thisrange >= minobrange) then azms: do j=1,strct_in_dbz(v,k)%num_beam !-Check to see if this is a missing observation nread=nread+1 if ( abs(strct_in_dbz(v,k)%field(i,j)) >= 99.0_r_kind ) then !--Extend no precip observations to missing data fields? ! May help suppress spurious convection if a problem. if (missing_to_nopcp) then strct_in_dbz(v,k)%field(i,j) = dbznoise num_m2nopcp = num_m2nopcp+1 else num_missing=num_missing+1 cycle azms !No reason to process the ob if it is missing end if end if if (l_limmax) then if ( strct_in_dbz(v,k)%field(i,j) > 60_r_kind ) then strct_in_dbz(v,k)%field(i,j) = 60_r_kind num_limmax=num_limmax+1 end if end if if (l_limmin) then if ( strct_in_dbz(v,k)%field(i,j) < 0_r_kind ) then strct_in_dbz(v,k)%field(i,j) = 0_r_kind num_limmin=num_limmin+1 end if end if !--Find observation height using method from read_l2bufr_mod.f90 this_stahgt=strct_in_dbz(v,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 !--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_dbz(v,k)%radlon clat0=cos(deg2rad*strct_in_dbz(v,k)%radlat) slat0=sin(deg2rad*strct_in_dbz(v,k)%radlat) thisazimuthr=(90.0_r_kind-strct_in_dbz(v,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 !-Check format of longitude and correct if necessary if(thislon>=r360) thislon=thislon-r360 if(thislon 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 netcdf use kinds, only: r_kind,r_double,i_kind,i_short,r_single use constants, only: zero,half,one,two,deg2rad,rearth,rad2deg, & one_tenth,r1000,r60,r60inv,r100,r400 use gridmod, only: tll2xy use obsmod, only: iadate use convinfo, only: nconvtype,ctwind,icuse,ioctype use mpimod, only: npe use read_l2bufr_mod, only : invtllv 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 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:: r360=360.0_r_kind integer(i_kind),parameter:: maxdat=17 ! Used in generating cdata array integer (i_kind):: iyear,imon,iday,ihour,imin,isec !--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,num_nopcp=0, & !counts numbadtime=0,num_badtilt=0, & num_badrange=0,num_m2nopcp=0, & num_noise=0,num_limmax=0 ,num_limmin=0 !--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,scount integer(i_kind),dimension(5) :: obdate real(r_kind) :: b,c,ha,epsh,h,aactual,a43,thistilt real(r_kind) :: thistiltr,selev0,celev0,thisrange,this_stahgt,thishgt real(r_kind) :: celev,selev,gamma,thisazimuthr,rlon0, & clat0,slat0,dlat,dlon,thiserr,thislon,thislat, & rlonloc,rlatloc,rlonglob,rlatglob,timeb,rad_per_meter real(r_kind) :: radartwindow real(r_kind) :: dbzerr,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_dbz(:,:) !---------SETTINGS FOR FUTURE NAMELIST---------! integer(i_kind) :: maxobrange=99900000 ! Range (m) *within* which to use observations - obs *outside* this range are not used integer(i_kind) :: minobrange=-999 ! Range (m) *outside* of which to use observatons - obs *inside* this range are not used real(r_kind) :: mintilt=0.0_r_kind ! Only use tilt(elevation) angles (deg) >= this number real(r_kind) :: maxtilt=20.0_r_kind ! Do no use tilt(elevation) angles (deg) >= this number logical :: missing_to_nopcp=.false. ! Set missing observations to 'no precipitation' observations -> dbznoise (See Aksoy et al. 2009, MWR) real(r_kind) :: dbznoise=2_r_kind ! dBZ obs must be >= dbznoise for assimilation logical :: l_limmax=.true. ! If true, observations > 60 dBZ are limited to be 60 dBZ. This is logical :: l_limmin=.true. ! If true, observations <0 dBZ are limited to be 0 dBZ. This is character (len=4) :: radarsite_nc character (len=256) vcpstr_nc !following the treatment on the precision issue for netcdf like in !wrf_netcdf_interface.F90 integer(i_kind) :: ncid,ierr,dimid1,dimid2,dimid3 integer(i_kind) :: varid1,varid2,varid3,varid4,varid6 integer(i_kind) :: pixel_x_varid,pixel_y_varid integer(i_kind) :: numazim_nc,numgate_nc,num_pixel_nc,real_num_pixel,vcp_nc real(r_single) :: elev_nc,firstgate_nc,lat_nc,lon_nc,height_nc integer(i_short),allocatable :: pixel_x_nc(:),pixel_y_nc(:) real(r_single), allocatable :: azimuth_nc(:),beamwidth_nc(:),azimspacing_nc(:),gatewidth_nc(:) real(r_single), allocatable :: nyquist_nc(:),obdata_nc(:,:),obdata_pixel_nc(:) real(r_single) nyquist_default_nc parameter(nyquist_default_nc=50.0_r_kind) logical l_pixel_unlimited integer(i_kind):: ipix integer(i_kind)::real_numpixel,start_nc(1),count_nc(1) !-Check if reflectivity is in the convinfo file and extract necessary attributes scount=0 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) dbzerr=5_r_kind !Ob error (dB) to use for radar reflectivity factor exit !Exit loop when finished with initial convinfo fields else if ( i==nconvtype ) then write(6,*) 'READ_dBZ: ERROR - OBSERVATION TYPE IS NOT PRESENT IN CONVINFO OR USE FLAG IS ZERO' write(6,*) 'READ_dBZ: ABORTTING read_dbz.f90 - NO REFLECTIVITY OBS READ!' return endif end do if (minobrange >= maxobrange) then write(6,*) 'MININMUM OB RANGE >= MAXIMUM OB RANGE FOR READING dBZ - PROGRAM STOPPING FROM READ_DBZ.F90' call stop2(400) end if !-next three values are dummy values for now nchanl=0 ilon=2 ilat=3 maxobs=2000000 !value taken from read_radar.f90 !--Allocate cdata_all array allocate(cdata_all(maxdat,maxobs)) lunrad=31 nvol=1 nelv=1 v=1;k=1 allocate(strct_in_dbz(nvol,nelv)) !!READ RADAR DATA ierr = NF90_OPEN(trim(infile),0,ncid) if (ierr /= nf90_noerr) call handle_err(ierr,"open") ierr = NF90_INQ_DIMID(ncid,'Azimuth',dimid1) if (ierr /= nf90_noerr) call handle_err(ierr,"Azimuth") ierr = NF90_INQ_DIMID(ncid,'Gate',dimid2) if (ierr /= nf90_noerr) call handle_err(ierr,"Gate") ierr = NF90_INQ_DIMID(ncid,'pixel',dimid3) if (ierr /= nf90_noerr) call handle_err(ierr,"Pixel number") ierr = NF90_INQ_VARID(ncid,'Azimuth',varid1) if (ierr /= nf90_noerr) call handle_err(ierr,"Azimuth") ierr = NF90_INQ_VARID(ncid,'BeamWidth',varid2) if (ierr /= nf90_noerr) call handle_err(ierr,"BeamWidth") ierr = NF90_INQ_VARID(ncid,'AzimuthalSpacing',varid3) if (ierr /= nf90_noerr) call handle_err(ierr,"azimuthalspacing") ierr = NF90_INQ_VARID(ncid,'GateWidth',varid4) if (ierr /= nf90_noerr) call handle_err(ierr,"gatewidth") ierr = NF90_INQ_VARID(ncid,'pixel_x',pixel_x_varid) ierr = NF90_INQ_VARID(ncid,'pixel_y',pixel_y_varid) ierr = NF90_INQ_VARID(ncid,'ReflectivityQC',varid6) if (ierr /= nf90_noerr) call handle_err(ierr,"ReflectivityQC") ierr = nf90_inquire_dimension(ncid, dimid1, len = numazim_nc) if (ierr /= nf90_noerr) call handle_err(ierr,"numazim data") ierr = nf90_inquire_dimension(ncid, dimid2, len = numgate_nc) if (ierr /= nf90_noerr) call handle_err(ierr,"numgate data") ierr = nf90_inquire_dimension(ncid, dimid3, len = num_pixel_nc) if (ierr /= nf90_noerr) call handle_err(ierr,"num_pixel_ncdata") if(num_pixel_nc<=0 ) then !unlimited size num_pixel_nc=numazim_nc*numgate_nc l_pixel_unlimited=.true. else real_num_pixel=num_pixel_nc l_pixel_unlimited=.false. endif ierr = NF90_GET_ATT(ncid,nf90_global,'Elevation',elev_nc) if (ierr /= nf90_noerr) call handle_err(ierr,"get elev") ierr = NF90_GET_ATT(ncid,nf90_global,'RangeToFirstGate',firstgate_nc) if (ierr /= nf90_noerr) call handle_err(ierr,"get firstgate") ierr = NF90_GET_ATT(ncid,nf90_global,'Latitude',lat_nc) if (ierr /= nf90_noerr) call handle_err(ierr,"get lat") ierr = NF90_GET_ATT(ncid,nf90_global,'Longitude',lon_nc) if (ierr /= nf90_noerr) call handle_err(ierr,"get lon") ierr = NF90_GET_ATT(ncid,nf90_global,'radarName-value',radarsite_nc) if (ierr /= nf90_noerr) call handle_err(ierr,"radarsite") ierr = NF90_GET_ATT(ncid,nf90_global,'vcp-value',vcpstr_nc) if (ierr /= nf90_noerr) call handle_err(ierr,"vcp") read(vcpstr_nc,*) vcp_nc ierr = NF90_GET_ATT(ncid,nf90_global,'Height',height_nc) if (ierr /= nf90_noerr) call handle_err(ierr,"height") !reverse order of dimensions as stated in ncdump: allocate(azimuth_nc(numazim_nc),beamwidth_nc(numazim_nc),azimspacing_nc(numazim_nc),gatewidth_nc(numazim_nc)) allocate(nyquist_nc(numazim_nc),obdata_nc(numgate_nc,numazim_nc)) allocate(obdata_pixel_nc(num_pixel_nc)) allocate(pixel_x_nc(num_pixel_nc)) allocate(pixel_y_nc(num_pixel_nc)) ierr = NF90_GET_VAR(ncid,varid1,azimuth_nc) if (ierr /= nf90_noerr) call handle_err(ierr,"azimuth data") ierr = NF90_GET_VAR(ncid,varid2,beamwidth_nc) if (ierr /= nf90_noerr) call handle_err(ierr,"beamwidth data") ierr = NF90_GET_VAR(ncid,varid3,azimspacing_nc) if (ierr /= nf90_noerr) call handle_err(ierr,"azimspacing data") ierr = NF90_GET_VAR(ncid,varid4,gatewidth_nc) if (ierr /= nf90_noerr) call handle_err(ierr,"gatewidth data") if(.not.l_pixel_unlimited) then ierr = NF90_GET_VAR(ncid,varid6,obdata_pixel_nc) real_numpixel=num_pixel_nc if (ierr /= nf90_noerr) call handle_err(ierr,"obdata_pixel data") ierr = NF90_GET_VAR(ncid,pixel_x_varid,pixel_x_nc) ierr = NF90_GET_VAR(ncid,pixel_y_varid,pixel_y_nc) else ierr=nf90_noerr ipix=1 start_nc=(/1/) count_nc=(/1/) ipix=1 do 255, while (ierr == nf90_noerr) start_nc(1)=ipix ierr = NF90_GET_VAR(ncid,varid6,obdata_pixel_nc(ipix:ipix),start=start_nc,count=count_nc) ierr = NF90_GET_VAR(ncid,pixel_x_varid,pixel_x_nc(ipix:ipix),start=start_nc,count=count_nc) ierr = NF90_GET_VAR(ncid,pixel_y_varid,pixel_y_nc(ipix:ipix),start=start_nc,count=count_nc) ipix=ipix+1 255 continue real_numpixel=ipix-2 endif ierr = NF90_CLOSE(ncid) if (ierr /= nf90_noerr) call handle_err(ierr,"close") do i=1,numazim_nc if ( (beamwidth_nc(i) /= beamwidth_nc(1)) .or. (gatewidth_nc(i) /= gatewidth_nc(1)) )then print *, "stopping: non-uniform scan" endif enddo read(infile(21:24),'(I4.4)')iyear read(infile(25:26),'(I2.2)')imon read(infile(27:28),'(I2.2)')iday read(infile(30:31),'(I2.2)')ihour read(infile(32:33),'(I2.2)')imin read(infile(34:35),'(I2.2)')isec do j=1,real_numpixel if(obdata_pixel_nc(j) < -999_r_kind) obdata_pixel_nc(j)=-999_r_kind enddo ! transform the read-in ob to the intermidate obs variables( radar obs to be used in GSI strct_in_dbz(v,k)%radid=radarsite_nc strct_in_dbz(v,k)%vcpnum=vcp_nc strct_in_dbz(v,k)%year=iyear ! to be defind from infile name strct_in_dbz(v,k)%month=imon strct_in_dbz(v,k)%day=iday strct_in_dbz(v,k)%hour=ihour strct_in_dbz(v,k)%minute=imin strct_in_dbz(v,k)%second=isec strct_in_dbz(v,k)%radlat=lat_nc strct_in_dbz(v,k)%radlon=lon_nc strct_in_dbz(v,k)%radhgt=height_nc strct_in_dbz(v,k)%fstgatdis =firstgate_nc strct_in_dbz(v,k)%gateWidth=gatewidth_nc(1) ! always the same ??) strct_in_dbz(v,k)%elev_angle=elev_nc strct_in_dbz(v,k)%num_beam=numazim_nc strct_in_dbz(v,k)%num_gate=numgate_nc na=strct_in_dbz(v,k)%num_beam nb=strct_in_dbz(v,k)%num_gate !******allocate arrays within radar data type**********! allocate(strct_in_dbz(v,k)%azim(na)) !******************************************************! strct_in_dbz(v,k)%azim(:)=azimuth_nc(:) ierror=0 fileopen: if (ierror == 0) then !Check to make sure file is open - will also fail if file does not exist. Closing endif at end of subroutine. !-Obtain analysis time in minutes since reference date call w3fs21(iadate,mins_an) !mins_an -integer number of mins snce 01/01/1978 ! w3movedat to help get a date from the time difference rmins_an=mins_an !convert to real number volumes: do v=1,nvol tilts: do k=1,nelv !--Check if observation fits within specified time window--! !-Find reference time of observation obdate(1)=strct_in_dbz(v,k)%year obdate(2)=strct_in_dbz(v,k)%month obdate(3)=strct_in_dbz(v,k)%day obdate(4)=strct_in_dbz(v,k)%hour obdate(5)=strct_in_dbz(v,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_dbz(v,k)%second*r60inv) !convert seconds to minutes and add to ob time !-Comparison is done in units of minutes timeb = rmins_ob-rmins_an ! now the window is controled by the preprocessing script starting from !5/14/2015 ! if(abs(timeb) > abs(radartwindow)) then ! numbadtime=numbadtime+1 ! cycle tilts !If not in time window, cycle the loop ! end if write(6,*) 'Processing obdate:',obdate,strct_in_dbz(v,k)%second if(abs(timeb) > 99999_r_kind) cycle !--Time window check complete--! thistilt=strct_in_dbz(v,k)%elev_angle if (thistilt <= maxtilt .and. thistilt >= mintilt) then pixel: do ipix=1,real_numpixel j=pixel_x_nc(ipix)+1 i=pixel_y_nc(ipix)+1 thisrange=strct_in_dbz(v,k)%fstgatdis + float(i-1)*strct_in_dbz(v,k)%gateWidth !-Check to make sure observations are within specified range if (thisrange <= maxobrange .and. thisrange >= minobrange) then nread=nread+1 if ( abs(obdata_pixel_nc(ipix)) >= 999.0_r_kind ) then !--Extend no precip observations to missing data fields? ! May help suppress spurious convection if a problem. if (missing_to_nopcp) then obdata_pixel_nc(ipix) = dbznoise num_m2nopcp = num_m2nopcp+1 else num_missing=num_missing+1 cycle pixel !No reason to process the ob if it is missing end if end if if (l_limmax) then if ( obdata_pixel_nc(ipix) > 60_r_kind ) then obdata_pixel_nc(ipix) = 60_r_kind num_limmax=num_limmax+1 end if end if if (l_limmin) then if ( obdata_pixel_nc(ipix) < 0_r_kind ) then obdata_pixel_nc(ipix) = 0_r_kind num_limmin=num_limmin+1 end if end if !-Special treatment for no-precip obs? !--Find observation height using method from read_l2bufr_mod.f90 this_stahgt=strct_in_dbz(v,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 !--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_dbz(v,k)%radlon clat0=cos(deg2rad*strct_in_dbz(v,k)%radlat) slat0=sin(deg2rad*strct_in_dbz(v,k)%radlat) thisazimuthr=(90.0_r_kind-strct_in_dbz(v,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 !-Check format of longitude and correct if necessary if(thislon>=r360) thislon=thislon-r360 if(thislon