! program sim_obs_conv ! ! Main program for creation of simulated conventional data for an OSSE ! ! Initial code by Ronald Errico January 2008 ! Assistence with bufr code by Meta Sienkiewicz ! Changes to print out by R. Errico August 2008 ! ! use m_kinds ! precision specification module for real variables use m_relhum ! transforms between q, r.h. and T_d use m_interp_NR ! NR read and interpolation module use m_bufr_rw ! Read/Write BUFR data ! ! link to libNCEP_bufr_r4i4.a library ! implicit none ! ! If ltest_all is true, then the program will use fake, low resolution NR ! data and fake observations and will print out many intermediate results ! for testing purposes. Alternatively, ltest=.true. or .false. can be ! specified at desired points in the code to print out only those portions. logical, parameter :: ltest_all=.false. ! ! Specify maximum dimensions for obs data arrays integer, parameter :: obs_ndim1=255 ! # of p-levels or channels integer, parameter :: obs_ndim2=2 ! # of fields of obs data (excl ps) integer, parameter :: obs_ndim3=1 ! # of p-lev info (excl obs values) integer, parameter :: obs_ndim4=8 ! # of scalar info for each obs ! obs_ndim1: maximum number of p-levels or channels for observations ! (currently > 255 not allowed sub. in read_write_obs) ! obs_ndim2: maximum number of observation fields in buffer (This only refers ! to obs derived from 3d-fields; e.g; 2 for u,v, but don't count ps) ! obs_ndim3: number of for each level or channel, besides the obs values ! themselves; e.g., p-levels, channels, and QC marks, etc. ! obs_ndim4: maximum number of single value report info (lat, lon etc.) ! integer, parameter :: n_errors=12 ! number of types of counters ! n_errors: number of types of counters for errors and observations processed ! logical :: ltest integer :: ierrors(n_errors) integer :: nlevs integer :: obs_nlevs integer :: obs_nfields integer :: nobs integer :: nfields2d ! number of required 2d fields integer :: nfields3d ! number of required 3d fields integer :: nfields ! sum of nfields2d and nfields3d character(len=4), allocatable :: f_names(:) ! list of field names ! real(rkind1) :: frac ! fraction of possible obs. values produced real(rkind3) :: obs_info(obs_ndim4) ! id,lon,lat,time,type,height,ps,... real(rkind3) :: obs_levels(obs_ndim1,obs_ndim3) ! plevs, QCflags, ... real(rkind3) :: obs_values(obs_ndim1,obs_ndim2) ! u,v or t,q, or channels ! ! Variables defined by Meta for modifying wind obs in prepbufr file integer :: luin, luout ! unit numbers integer :: argc integer(4) :: iargc integer :: ireadmg integer :: ireadsb integer :: idate ! synoptic date/time YYYYMMDDHH character(len=5) :: d_type ! either MASS or WIND character(len=120) :: inputfile character(len=120) :: outputfile character(len=8) :: subset ! name of current BUFR subset character(len=8) :: psubset ! name of previous BUFR subset logical :: check_type logical :: modtype ! flag indicats whether type of report to modify real(rkind3) :: bmiss ! value of flag to indicate missing value logical, parameter :: lsim=.true. ! program simulates observations data luin /8/, luout /9/ ! do not use units 10-12 ! ! Read and check arguments argc = iargc() if (argc .lt. 3) then print *,'usage must be: cdata_Main.x d_type inputbufr outputbufr' stop endif call GetArg( 1_4, d_type) call GetArg( 2_4, inputfile) call GetArg( 3_4, outputfile) print *,' ' if ( (d_type == 'WIND_') .or. (d_type == 'MASS_') ) then print *,'Begin processing for type=',d_type else print *,'Specification of d_type=',d_type,' neither ', & 'WIND_ nor MASS_ as required' stop endif ! ! Specify NR fields to be interpolated if (d_type=='WIND_') then nfields2d=1 nfields3d=2 nfields=nfields2d+nfields3d allocate (f_names(nfields)) f_names(1)='pres' f_names(2)='uwnd' f_names(3)='vwnd' elseif (d_type=='MASS_') then nfields2d=2 nfields3d=2 nfields=nfields2d+nfields3d allocate (f_names(nfields)) f_names(1)='pres' f_names(2)='zsfc' f_names(3)='temp' f_names(4)='sphu' endif ! ! Rules for setting names for setting f_names to set of NR fields to be used: ! a) choose field names from list appearing at beginning of module m_interp_NR ! These may be in any order, except for the following rules: ! a) all 2D fields should preceed 3D fields ! b) 'uwnd' should be followed by 'vwnd' ! c) 'sphu' should be preceeded 'temp' ! ! Set or initialize some variables ierrors(:)=0 bmiss = 10.e10 ltest=ltest_all ! ! Set dimensions of grid and allocate module arrays call Setup_m_interp (nfields2d,nfields3d,nlevs,ltest) ! ! Read in and set up either NR or (if testing) fake data call Setup_NR_fields (f_names,ltest) ! ! Use fake data if testing if (ltest) then do nobs=1,4 call read_obs_test (obs_ndim1,obs_ndim2,obs_ndim3, & obs_ndim4,obs_nlevs,obs_nfields, & nobs,obs_info,obs_levels,obs_values) call Interp_NR_fields (obs_ndim1,obs_ndim2,obs_ndim3, & obs_ndim4,obs_nlevs,obs_nfields, & nobs,n_errors,ierrors,bmiss, & obs_info,obs_levels,obs_values,ltest) enddo endif ! Setup files for reading and writing observations if (.not.ltest) then ! open(unit=luin,file =trim(inputfile), form='unformatted') print *,' ' print ('(3a,i3)'),' input file=',trim(inputfile), & ' opened on unit=',luin open(unit=luout,file=trim(outputfile),form='unformatted') print *,' ' print ('(3a,i3)'),' output file=',trim(outputfile), & ' opened on unit=',luout call openbf(luin, 'IN ',luin) call openbf(luout,'OUT',luin) ! psubset = '' modtype = .false. ! ! Begining of loop to read observations on existing file nobs=0 call datelen (10) ! returns a 10 digit value of idate YYYYMMDDHH do while (ireadmg(luin,subset,idate).eq. 0) ! If new data source (new 'subset') if (subset .ne. psubset ) then print ('(3a,i10)'),'Processing subset ',subset, & ' for datetime ',idate psubset = subset endif ! set flag 'modtype' if we want to modify this data type (i.e. if it has 't') modtype = check_type (subset, d_type) if ( modtype ) then call openmb(luout,subset,idate) ! do while ( ireadsb(luin) .eq. 0 ) ! ! End of begin loop instructions from Meta if (d_type == 'WIND_' ) then call read_write_obs_uv (obs_ndim1,obs_ndim2,obs_ndim3, & obs_ndim4,obs_nlevs,obs_nfields,nobs, & ierrors(9),luin,luout,ltest,lsim,bmiss, & obs_info,obs_levels,obs_values,.true.) else call read_write_obs_tq (obs_ndim1,obs_ndim2,obs_ndim3, & obs_ndim4,obs_nlevs,obs_nfields,nobs, & ierrors(9),luin,luout,ltest,lsim,bmiss, & obs_info,obs_levels,obs_values,.true.) endif ! ! Create new observations by interpolation from NR call interp_NR_fields (obs_ndim1,obs_ndim2,obs_ndim3, & obs_ndim4,obs_nlevs,obs_nfields, & nobs,n_errors,ierrors,bmiss, & obs_info,obs_levels,obs_values,.false.) ! ! Write out observation buffer if (d_type == 'WIND_' ) then call read_write_obs_uv (obs_ndim1,obs_ndim2,obs_ndim3, & obs_ndim4,obs_nlevs,obs_nfields,nobs, & ierrors(9),luin,luout,ltest,lsim,bmiss, & obs_info,obs_levels,obs_values,.false.) else call read_write_obs_tq (obs_ndim1,obs_ndim2,obs_ndim3, & obs_ndim4,obs_nlevs,obs_nfields,nobs, & ierrors(9),luin,luout,ltest,lsim,bmiss, & obs_info,obs_levels,obs_values,.false.) endif ! ! End loop over buffer records enddo ! loop over do while ( ireadsb(luin) .eq. 0 ) call closmg(luout) print *,'QQQQ',obs_nlevs,obs_nfields do argc=1,obs_nlevs print ('(i2,4e15.6)'),argc, & obs_levels(argc,1),obs_values(argc,1:obs_nfields) enddo endif ! check on modtype enddo ! loop over (ireadmg(luin,subset,idate).eq. 0) call closbf(luout) ! endif ! test of ltest ! ! Print out summary table of number of obs reports processed print *,' ' print *,' SUMMARY TABLE:' print ('(i8,2a)'),nobs,' observation reports read' print ('(i8,a)'),ierrors(10), & ' number of reports without data or not requested data types' print ('(i8,a)'),ierrors(11), & ' number of reports having some data of requested types' frac=real(ierrors(11))/real(ierrors(10)+ierrors(11)) print ('(f8.5,a)'),frac, & ' fraction of reports read having requested data types' print ('(i8,a)'),ierrors(8),' number of observation values considered' frac=1.-real(ierrors(7))/real(ierrors(8)) print ('(f8.5,a)'),frac, & ' fraction of obs values simulated vs. read for requested types' ! ! Print out error count summary ! ierrors(6)=nobs is number of observation reports read in ierrors(n_errors)=ierrors(1)+ierrors(2)+ierrors(3)+ierrors(4)+ & ierrors(5)+ierrors(6)+ierrors(7)+ierrors(9)+ & ierrors(10) if (ierrors(n_errors) > 0) then print *,' ' print *,' Summary of bad observations or other errors detected:' print ('(i8,a)'),ierrors(1),' observation reports found where ttmax' print ('(i8,a)'),ierrors(3), & ' observation reports found where longitude out of range' print ('(i8,a)'),ierrors(4), & ' observation reports found where latitude out of range' print ('(i8,a)'),ierrors(5), & ' observation values found where obs_plev < ptop' print ('(i8,a)'),ierrors(6), & ' observation values where obs_plev > ps lowest level' print ('(i8,a)'),ierrors(7), & ' observation values ignored for various detected problems' print ('(i8,a)'),ierrors(9), & ' errors detected in writing buffer records' endif ! deallocate (f_names) call clean_m_interp ! deallocate grid and field arrays ! print *,' ' print *,'Program completed' ! end program sim_obs_conv !