! program sim_obs_rad ! ! Main program for creation of simulated radiance data for an OSSE ! ! Initial code by Ronald Errico January 2008 ! Assistence with bufr code by Meta Sienkiewicz and Hui-Chun Liu. ! Interfacing with CRTM by Runhua Yang January/Feb 2008 ! Minor corrections and changes to print out by R. Errico August 2008 ! Additions to simulate MSU and to include cloud liquid water Errico Sept 2008 ! ! ! program history log: ! 2008-01-xx errico - initial code ! 2008-02-xx yang - add crtm interface ! 2008-09-17 errico - add cloud liquic water and msu ! 2008-09-19 yang - modify to include msu. Test ! 2008-10-27 errico - change loop structure to enable easy mpi 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 use m_rdata_boxes ! thinning boxes use m_cloud ! determines if radiances are cloud or precip affected use crtm_module ! only the channel info use m_interface_crtm ! interface module to call crtm_forwardmodel ! ! 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. ! ltest_sample .true. cause print of some output for comparing results logical, parameter :: ltest_all=.false. logical, parameter :: ltest_sample=.true. ! ! Specify maximum dimensions for obs data arrays integer, parameter :: obs_ndim1=350 integer, parameter :: obs_ndim2=2 integer, parameter :: obs_ndim3=2 integer, parameter :: obs_ndim4=3*12+25 ! number of AQUA/AIRS info ! obs_ndim1: maximum number of p-levels or channels for observations ! obs_ndim2: maximum number of observation data info in buffer ! (2 includes both brightness temperature and quality mark) ! obs_ndim3: number of channel info ! (2 includes both channel number and log reciprical wavelength) ! obs_ndim4: maximum number of single value report info (lat, lon etc.) ! integer, parameter :: n_errors=12 ! number of types of counters integer, parameter :: n_clouds=6 ! number of counted cloud levels integer, parameter :: n_sfc_types=3 ! number of sfc types distinguished ! logical :: l_box ! true if info saved in thinning box logical :: lerror1, lerror2 ! true if errors detected by subroutine logical :: lfirst logical :: lmsg1 ! .true if first message on file has 0 recs logical :: ltest ! true if some info to be printed for test integer :: box_index, box_index2, box_s2 integer, allocatable :: box_nums(:,:) ! counters of non-empty boxes integer, allocatable :: box_skip(:) ! order of types having some data integer :: i, i1, ix1, ix2, ix3 integer*8 :: idatetime integer :: ierrors(n_errors) integer :: ierrors2(n_errors) integer :: ierrors3(n_errors) integer :: i_obs integer :: j_obs integer :: j_obs_M integer :: j_obs_last integer :: j_obs_span integer :: lcldtop integer :: nlevs integer :: nboxes ! number of thinning boxes integer :: n_channels_crtm integer :: n_mesg(2) ! counter of messages in files integer :: n_mesg_next ! counter of next message to be read integer :: n_subtypes ! number of subset sat or instr ID #s integer :: obs_n_channels integer :: obs_n_data integer :: nfields2d ! number of required 2d fields integer :: nfields3d ! number of required 3d fields integer :: nfields ! sum of nfields2d and nfields3d integer :: nobs ! counter for number of observations read integer :: nobs2 integer :: n_processors ! number of processors used (for MPI) integer :: n_spot ! size of observation instrument report header integer :: n_spot2 ! size of 'AQUASPOT' report header integer :: prof_ndim4 integer :: n_subtype ! max number of satellite subtypes on file integer :: sfc_type ! flag for 1=ocean, 2=ice, 3=land integer :: sfc_types(n_sfc_types,2) ! counter for profile sfc types integer :: subtype ! ID # for either satellite or intrument type ! character(len=4), allocatable :: f_names(:) ! list of field names ! real(rkind1) :: box_size ! thinning box size (km) real(rkind1) :: frac ! fraction of possible obs. values produced real(rkind1) :: cloudy(n_clouds,3) ! for counting "cloud" effects real(rkind1) :: pcldtop real(rkind1) :: psfc real(rkind3) :: obs_info(obs_ndim4) ! id,lon,lat,time,type,scan, etc. real(rkind3) :: obs_channels(obs_ndim1,obs_ndim3) ! channel id and freq real(rkind3) :: obs_values(obs_ndim1,obs_ndim2) ! B-temp and qc flag real(rkind1) :: soza ! solar zenith angle ! ! prof_info contains header for a single satellite instrument augmented by ! 4 extra values: (1) either counter of numbers of obs that fell in the ! thinning box containing this obs or, in the case of AMSUA/AQUA, a flag ! indicating whether the header has missing data or not; (2) latitude; ! (3) longitude; (4) hours before or after synoptic time of file. real(rkind1), allocatable :: prof_info(:) real(rkind1), allocatable :: prof_2d(:) real(rkind1), allocatable :: prof_3d(:,:) real(rkind1), allocatable :: prof_plevs(:) ! ! A saved copy of a profile in case one for AMSUA/AQUA is invalid real(rkind1), allocatable :: prof_info_S(:) real(rkind1), allocatable :: prof_2d_S(:) real(rkind1), allocatable :: prof_3d_S(:,:) real(rkind1), allocatable :: prof_plevs_S(:) ! ! Arrays for saving information for multiple profiles for MPI logical, allocatable :: lerror2_M(:) integer, allocatable :: lcldtop_M(:) integer, allocatable :: sfc_type_M(:) real(rkind1), allocatable :: prof_info_M(:,:) real(rkind1), allocatable :: prof_plevs_M(:,:) real(rkind1), allocatable :: prof_2d_M(:,:) real(rkind1), allocatable :: prof_3d_M(:,:,:) real(rkind1), allocatable :: pcldtop_M(:) real(rkind1), allocatable :: soza_M(:) real(rkind3), allocatable :: crtmBT_M(:,:) real(rkind3), allocatable :: obs_info_M(:,:) real(rkind3), allocatable :: obs_values_M(:,:,:) ! integer :: luin, luout, luin_tab, lrc ! unit numbers integer :: i_return integer :: argc integer(4) :: iargc integer :: ireadmg integer :: ireadsb integer :: idate ! date/time YYYYMMDDHH read on file integer :: idate1 ! copy of time on first message character(len=5) :: d_type ! data type character(len=120) :: c_datetime character(len=120) :: rcfile character(len=120) :: inputfile character(len=120) :: outputfile character(len=120) :: output_name character(len=8) :: subset ! name of current BUFR subset character(len=8) :: subset1 ! copy of first subset on file 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 ! type(crtm_channelinfo_type):: channelinfo integer :: nchan integer :: error_status ! data luin /8/, luout /9/, luin_tab /15/ ! do not use units 10-12 data lrc /16/ ! do not use units 10-12 ! ! print *,' ' print *,'BEGIN PROGRAM sim_obs_rad' ! ! Read and check arguments argc = iargc() if (argc .lt. 5) then print *,' usage must be: sim_rad.x d_type datetime', & ' rcfile inputbufr outputbufr' stop endif call GetArg( 1_4, d_type) call GetArg( 2_4, c_datetime) call GetArg( 3_4, rcfile) call GetArg( 4_4, inputfile) call GetArg( 5_4, outputfile) print *,' ' if ( (d_type == 'AIRS_') .or. (d_type == 'HIRS2') .or. & (d_type == 'HIRS3') .or. (d_type == 'AMSUA') .or. & (d_type == 'AMSUB') .or. (d_type == '_MSU_')) then print *,'Begin processing for type=',d_type else print *,'Specification of d_type=',d_type,' neither ', & 'AIRS_, HIRS2, HIRS3, AMSUA, AMSUB, nor _MSU_ as required' stop endif read (c_datetime(1:10),'(i10)') idatetime ! ! Specify NR fields to be interpolated ! The order of f_names must be consistent with prof_2d and prof_3d indexes ! in this program appearing below. if (d_type(1:4)=='HIRS') then ! fields for IR calculations nfields2d=9 nfields3d=3 nfields=nfields2d+nfields3d allocate (f_names(nfields)) f_names(1)='pres' f_names(2)='zsfc' f_names(3)='hcld' f_names(4)='mcld' f_names(5)='lcld' f_names(6)='u10m' f_names(7)='v10m' f_names(8)='ismk' f_names(9)='almk' f_names(10)='temp' f_names(11)='sphu' f_names(12)='ozon' output_name=outputfile ! elseif (d_type(1:4)=='AIRS') then ! for both IR and MW calculations nfields2d=11 nfields3d=3 nfields=nfields2d+nfields3d allocate (f_names(nfields)) f_names(1)='pres' f_names(2)='zsfc' f_names(3)='hcld' f_names(4)='mcld' f_names(5)='lcld' f_names(6)='u10m' f_names(7)='v10m' f_names(8)='ismk' f_names(9)='almk' f_names(10)='conp' f_names(11)='rain' f_names(12)='temp' f_names(13)='sphu' f_names(14)='ozon' output_name=trim(outputfile)//'_temp' ! elseif (d_type(2:4)=='MSU') then ! fields for MW calculations nfields2d=9 nfields3d=3 nfields=nfields2d+nfields3d allocate (f_names(nfields)) f_names(1)='pres' f_names(2)='zsfc' f_names(3)='sktp' f_names(4)='u10m' f_names(5)='v10m' f_names(6)='ismk' f_names(7)='almk' f_names(8)='conp' f_names(9)='rain' f_names(10)='temp' f_names(11)='sphu' f_names(12)='ozon' output_name=outputfile else print ('(3a)'),' Option d_type=',d_type,' not supported' stop 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 n_mesg(:)=0 lmsg1=.false. ierrors(:)=0 ierrors2(:)=0 bmiss = 10.e10 ltest=ltest_all if (d_type=='AIRS_') then n_spot=12 ! # values in instrument header n_spot2=25 ! # values in "AQUASPOT" header n_subtypes=3 ! 3 instruments on AIRS/AQUA file else n_spot=16 ! # values in spot header n_spot2=0 ! not used n_subtypes=5 ! max of 5 satellites for AMSU-B endif ! ! 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) ! ! Read in resource file for specifying cloud, precip., or surface effects call Set_cloud (idatetime,nfields2d,nlevs,lrc,rcfile,d_type, & f_names,box_size) ! ! Initialize array for counting the numbers of observations affected by ! clouds, precip., or the surface. The effects are distinguished by counting ! the numbers of profiles with effective radiating surface elevated at a ! level within some range of values of sigma=p/ps. In the case here, there ! are n_clouds-1 ranges, each of width delta sigma = 1/(n_cloud-1) cloudy(:,2:3)=zero_k1 ! initialize counter do i=1,n_clouds cloudy(i,1)=real(n_clouds-i)/real(n_clouds-1) ! edges of ranges enddo sfc_types(:,:)=0 ! initialize counter ! ! Allocate some arrays prof_ndim4=n_spot+4 allocate (prof_info(prof_ndim4)) allocate (prof_2d(nfields2d)) allocate (prof_3d(nlevs,nfields3d)) allocate (prof_plevs(nlevs)) ! ! Set thinning box geometry and arrays call set_rdata_boxes (box_size,n_spot,prof_ndim4,nboxes,n_subtypes, & n_spot2) ! ! Open AIRS bufr table if required if (d_type=='AIRS_') then open(unit=luin_tab,file ='airs_bufr_table', form='formatted') print *,' ' print ('(3a,i4)'),' input file=','airs_bufr_table', & ' opened on unit=',luin_tab else luin_tab=luin ! assumes bufr table appended to input file endif ! ! Begin first pass through data. This pass essentially reads all ! observation reports and retains a thinned subset of header information: ! (1) read each report ! (2) interpolate 2-d information from the NR to the report location ! (3) use the 2-d information to assign cloud or surface effects ! (4) save the best information falling into each thinning box ! ! Open input bufr data file open(unit=luin,file =trim(inputfile), form='unformatted') print *,' ' print ('(3a,i3)'),' input file=',trim(inputfile), & ' opened on unit=',luin call openbf(luin, 'IN ',luin_tab) ! ! Begin loop to read observations from existing file nobs=0 psubset = '' 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 date ',idate psubset = subset if (nobs==0) then ! Save first value of message header with info idate1=idate subset1=subset endif endif ! modtype = check_type (subset, d_type) if (modtype) then do while ( ireadsb(luin) .eq. 0 ) ! if (d_type == 'AIRS_' ) then call read_write_obs_airs (obs_ndim1,obs_ndim2,obs_ndim3, & obs_ndim4,obs_n_channels,obs_n_data,nobs, & ierrors(9),luin,luout,ltest,lsim,bmiss, & obs_info,obs_channels,obs_values,.true.) else call read_write_obs_amsu_hirs (obs_ndim1,obs_ndim2,obs_ndim3, & obs_ndim4,obs_n_channels,obs_n_data,nobs, & ierrors(9),luin,luout,ltest,lsim,bmiss, & obs_info,obs_channels,obs_values,.true.) endif ! ! Test some obs header info, compute departure from snoptic time, ! and fill the array prof_info call fill_prof_info (idatetime,n_errors,ierrors, & obs_info(1:n_spot),d_type,prof_info,lerror1) call find_rad_subset (prof_ndim4,prof_info,d_type, & ierrors(5),lerror1,box_index2) ! ! If AIRS, check that report also contains a valid AMSUA header ! (actually, only check second value of header) if (d_type == 'AIRS_') then if (obs_info(n_spot+2) > 1.e8) then lerror1=.true. endif endif ! ! Interpolate 2-d NR fields used to assign observation quality ! (e.g., cloud or surface effects) if (.not.lerror1) then ! continue processing call find_box_index (prof_info(2),prof_info(3),box_index) call interp_fields_for_rad (nlevs,nfields3d,nfields2d, & prof_ndim4,n_spot,nobs,n_errors,ierrors, & prof_info,prof_2d,prof_3d,prof_plevs, & .false.,lerror2,ltest) else ! do not process; just flag as error lerror2=.true. endif ! test on .not.lerror1 ! ! Determine observation quality based on cloud/surface/precip conditions ! Then save the best quality information found for each thinning box if (.not.lerror2) then call cloud (nlevs,nfields2d,prof_2d,prof_plevs, & lcldtop,pcldtop) call put_box_info (box_index,box_index2,lcldtop,pcldtop, & prof_info,l_box) ! ! Save AMSU and HSB header embedded with AIRS information in report ! Use same values of lcldtop and pcldtop as for AIRS here, so that ! same report is retained if (l_box .and. (d_type == 'AIRS_')) then if (n_spot2 > 0) then call put_box_spot_info (box_index, & obs_info(3*n_spot+1:3*n_spot+n_spot2)) endif call fill_prof_info (idatetime,n_errors,ierrors2, & obs_info(n_spot+1:2*n_spot),d_type,prof_info,lerror1) call put_box_info (box_index,2,lcldtop,pcldtop, & prof_info,l_box) call fill_prof_info (idatetime,n_errors,ierrors2, & obs_info(2*n_spot+1:3*n_spot),d_type,prof_info,lerror1) call put_box_info (box_index,3,lcldtop,pcldtop, & prof_info,l_box) endif ! endif ! test on lerror2 ! ! End loops to read in observation buffers enddo ! loop over do while ( ireadsb(luin) .eq. 0 ) if (nobs==0) lmsg1=.true. ! indicates first message has 0 reports endif ! check on modtype n_mesg(1)=n_mesg(1)+1 ! count messages read in original file enddo ! loop over (ireadmg(luin,subset,idate).eq. 0) ierrors(6)=nobs ! ! Determine the numbers of observations saved for each subtype ! And determine the sequence of indexes for box subtypes that have info. ! box_skip(i)=0 if no boxes filled for that type ! box_nums(:,2) are numbers of filled boxes for each type ! box_nums(:,1) is counter of profiles processed for each type ! (the latter is used to select sample of plots to print out for testing) allocate (box_skip(n_subtypes)) allocate (box_nums(n_subtypes,2)) call set_box_skip (nboxes,n_subtypes,box_skip,box_nums(:,2)) box_nums(:,1)=0 ! ! End first pass through data ! ! Begin second pass through data. This pass essentially reads one ! observation report for each thinning box and then: ! (1) interpolates the 3D fields to create an atmospheric column profile ! (2) calls the CRTM ! (3) writes out the simulated observation report ! ! If this is for AIRS data, results are written to a temporary file, whose ! data are augmented later when the AMSU-A AQUA data are added. At that ! point, the requested output file is produced. ! ! Open output bufr data file open(unit=luout,file=trim(output_name),form='unformatted') print *,' ' print ('(3a,i3)'),' output file=',trim(output_name), & ' opened on unit=',luout call openbf(luout,'OUT',luin_tab) ! ! Compress messages in all BUFR files written (so several records fit in 1 msg) call cmpmsg('Y') ! ! Set some variables to start second loop call setupf_names (nfields2d,nfields3d,nfields,f_names) if (d_type == 'AIRS_') then n_channels_crtm=281 elseif (d_type (1:4) == 'HIRS' ) then n_channels_crtm=19 else n_channels_crtm=obs_n_channels print *, 'd_type=', d_type, n_channels_crtm endif nchan=n_channels_crtm ! ! For AIRS processing, need to save some profile information for AMSUA if (d_type == 'AIRS_') then allocate (prof_info_S(prof_ndim4)) allocate (prof_2d_S(nfields2d)) allocate (prof_3d_S(nlevs,nfields3d)) allocate (prof_plevs_S(nlevs)) endif ! ! Allocate arrays for communicating channel information with CRTM !MMMMMM (1) set n_processors to # of processors to be used in expensive loop n_processors=10 allocate (prof_info_M(prof_ndim4,n_processors)) allocate (prof_plevs_M(nlevs,n_processors)) allocate (prof_2d_M(nfields2d,n_processors)) allocate (prof_3d_M(nlevs,nfields3d,n_processors)) allocate (lcldtop_M(n_processors)) allocate (pcldtop_M(n_processors)) allocate (soza_M(n_processors)) allocate (crtmBT_M(n_channels_crtm,n_processors)) allocate (sfc_type_M(n_processors)) ! print *,' ' print ('(3a,i10)'),' 2nd loop Processing subset ',subset1, & ' for date ',idate1 ! ! If the first message in the original file had only the date/time stamp with ! no reports, then in the file to be created, do likewise so that GSI reads ! the file properly. call openmg(luout,subset1,idate1) ! add subset and idate to new file call minimg(luout,0) ! add minutes=00 to idate1 write if (lmsg1) call closmg(luout) ! ! Begin loop over non-empty thinning boxes do box_s2=1, n_subtypes ! loop over possible subtypes box_index2=box_skip(box_s2) ! index for current present subtype ! ! Exit if no more subtypes to be processed for this pass of data if ((box_index2 == 0) .or. & ((d_type == 'AIRS_') .and. (box_s2 == 2)) ) exit ! ! Loop over all thinned observations of current subtype in sets of ! n_processors each (except for the last set, that may be fewer) do i_obs=1, box_nums(box_index2,2), n_processors j_obs_last=min (i_obs + n_processors - 1, box_nums(box_index2,2)) j_obs_span=j_obs_last - i_obs + 1 ! number of obs current set ! ! Loop to fill arrays holding sets of profiles to pass to CRTM do j_obs=i_obs, j_obs_last ! ! Retrieve stored information from next thinning box call get_box_info (j_obs,box_index2,prof_info) lcldtop =nint(prof_info(prof_ndim4-2)) pcldtop =prof_info(prof_ndim4-1) ! ! Interpolate nature-run fields create profiles at obs. locations call interp_fields_for_rad (nlevs,nfields3d,nfields2d, & prof_ndim4,n_spot,j_obs,n_errors,ierrors, & prof_info,prof_2d,prof_3d,prof_plevs, & .true.,lerror2,ltest) if (lerror2) then print *,' ' print *,'Problem when interpolating during 2nd loop' print *,'This should not be possible since data already',& ' checked when thinning' stop endif ! ! Accumulate counts for each cloud range if (lcldtop == nlevs) then cloudy(1,2)=cloudy(1,2)+one_k1 ! cloud free observation else do i=2,n_clouds if ((pcldtop/prof_2d(1) < cloudy(i-1,1)) & .and. (pcldtop/prof_2d(1) >= cloudy(i,1))) then cloudy(i,2)=cloudy(i,2)+one_k1 endif enddo endif ! ! Get solar zenith angle. For AIRS it is in the 'AQUASPOT' header, ! but for the NOAA satellites it is in the instrument headers. if (d_type == 'AIRS_' ) then call get_box_spot_info (j_obs, & obs_info(3*n_spot+1:3*n_spot+n_spot2)) soza=obs_info(3*n_spot+6) else soza=prof_info(12) endif ! ! Copy profile of each obs in set to the array of multiple profiles j_obs_M=j_obs-i_obs+1 prof_info_M(:,j_obs_M)=prof_info(:) prof_plevs_M(:,j_obs_M)=prof_plevs(:) prof_2d_M(:,j_obs_M)=prof_2d(:) prof_3d_M(:,:,j_obs_M)=prof_3d(:,:) lcldtop_M(j_obs_M)=lcldtop pcldtop_M(j_obs_M)=pcldtop soza_M(j_obs_M)=soza ! enddo ! end first loop over j_obs ! ! If this is first occurance of subtype, initialize some values if (i_obs == 1) then ! if (d_type == 'AIRS_') then subtype=nint(prof_info(1)) ! instrument ID # else subtype=nint(prof_info(14)) ! satellite ID # endif lfirst=.true. ! ! Initialize multiple copies of the CRTM !MMMM (2) MPI this loop do j_obs_M=1,1 !MMMM in MPI version, must be =1,j_obs_span error_status = crtm_destroy(channelinfo) if (error_status /= success) then print *,'SIM_OBS_RAD: ***ERROR*** crtm_destroy ', & 'error_status=',error_status endif call InitGetbt_crtm (prof_ndim4,nlevs,nfields2d,nfields3d, & prof_info_M(:,j_obs_M),prof_plevs_M(:,j_obs_M), & prof_2d_M(:,j_obs_M),prof_3d_M(:,:,j_obs_M), & lcldtop_M(j_obs_M),pcldtop_M(j_obs_M), & lfirst,d_type,subtype,crtmBT_M(:,j_obs_M),nchan, & channelinfo,soza_M(j_obs_M),sfc_type_M(j_obs_M)) enddo ! print *, 'Old channelinfo destroyed: new subtype= ', subtype endif ! test on i_obs == 1 ! ! Call CRTM to create sets of observations from sets of profiles lfirst=.false. !MMMM (3) MPI this loop do j_obs_M=1,j_obs_span call InitGetbt_crtm (prof_ndim4,nlevs,nfields2d,nfields3d, & prof_info_M(:,j_obs_M),prof_plevs_M(:,j_obs_M), & prof_2d_M(:,j_obs_M),prof_3d_M(:,:,j_obs_M), & lcldtop_M(j_obs_M),pcldtop_M(j_obs_M), & lfirst,d_type,subtype,crtmBT_M(:,j_obs_M),nchan, & channelinfo,soza_M(j_obs_M),sfc_type_M(j_obs_M)) enddo ! ! Write report of simulated obs do j_obs=i_obs, j_obs_last ! ! Copy simulated for obs values at one location from set of multiple obs j_obs_M=j_obs - i_obs + 1 obs_values(1:n_channels_crtm,1)= & crtmBT_M(1:n_channels_crtm,j_obs_M) sfc_type=sfc_type_M(j_obs_M) sfc_types(sfc_type,1)=sfc_types(sfc_type,1)+1 prof_info(:)=prof_info_M(:,j_obs_M) obs_info(1:n_spot)=prof_info_M(1:n_spot,j_obs_M) ierrors(7)=ierrors(7)+1 ! number of error-free thinned obs. ! call openmb(luout,subset1,idate1) ! ! Special instructions for AIRS because report also contains other instruments if (d_type == 'AIRS_' ) then ! ! Save a single profile for substitute of any possible invalid ones ! encountered in data pass 3. if (ierrors(7)==1) then prof_info_S(:)=prof_info(:) prof_2d_S(:)=prof_2d(:) prof_3d_S(:,:)=prof_3d(:,:) prof_plevs_S(:)=prof_plevs(:) endif ! ! Save headers for AMSUA, HSB, and AQUA-SPOT call get_box_spot_info (j_obs, & obs_info(3*n_spot+1:3*n_spot+n_spot2)) call get_box_info (j_obs,2,prof_info) obs_info(n_spot+1:2*n_spot)=prof_info(1:n_spot) call get_box_info (j_obs,3,prof_info) obs_info(2*n_spot+1:3*n_spot)=prof_info(1:n_spot) call read_write_obs_airs (obs_ndim1,obs_ndim2,obs_ndim3, & obs_ndim4,obs_n_channels,obs_n_data,ierrors(7), & ierrors(9),luin,luout,ltest,lsim,bmiss, & obs_info,obs_channels,obs_values,.false.) ! else ! not an AIRS report call read_write_obs_amsu_hirs (obs_ndim1,obs_ndim2, & obs_ndim3, & obs_ndim4,obs_n_channels,obs_n_data,ierrors(7), & ierrors(9),luin,luout,ltest,lsim,bmiss, & obs_info,obs_channels,obs_values,.false.) endif ! ! Print sample of data if testing code box_nums(box_index2,1)=box_nums(box_index2,1)+1 if (ltest_sample .and. & mod(box_nums(box_index2,1),box_nums(box_index2,2)/4)==1) then if (n_spot2 == 0) then ix1=obs_ndim4 ! no "spot" header present ix2=ix1 else ix1=n_spot*n_subtypes+1 ! AIRS "spot" header present ix2=ix1+n_spot2-1 endif ix3=ix2-ix1+1 if (d_type == 'AIRS_') then ! re-copy AIRS header prof_info(:)=prof_info_M(:,j_obs_M) endif call print_sample (d_type,box_nums(box_index2,1),subtype, & ierrors(7),prof_ndim4,n_spot,n_spot2, & ix3,n_channels_crtm,sfc_type,prof_info, & obs_info(1:n_spot),obs_info(ix1:ix2), & obs_values(1:n_channels_crtm,1)) endif ! test on ltest_sample, etc. enddo ! second loop over j_obs ! ! End loops to create observations enddo ! loop over i_obs (all obs for 1 subtype) ! Accumulate total number of empty thinning boxes ierrors(8)=ierrors(8)+nboxes-box_nums(box_index2,1) enddo ! loop over box_s2 (all possible subtypes) call closmg(luout) ! final close out of messages ! ! End loop over buffer records call closbf(luin) call closbf(luout) call clean_rdata_boxes ! deallocate thinning boxes call clean_cloud ! deallocate table of thinning criteria deallocate (crtmBT_M) if ((d_type /= 'AIRS_') .and. (luin /= luin_tab)) call closbf(luin_tab) ! ! ! End second pass through data ! ! Check for need of 3rd pass (required for AMSUA/AQUA only) ! if (d_type == 'AIRS_') then ! ! Re-open temporary AIRS output file open(unit=luin,file=trim(output_name),form='unformatted') print *,' ' print ('(3a,i3)'),' temporary file=',trim(output_name), & ' re-opened for input on unit=',luin call openbf(luin,'IN',luin_tab) ! ! Open actual output bufr data file open(unit=luout,file=trim(outputfile),form='unformatted') print *,' ' print ('(3a,i3)'),' output file=',trim(outputfile), & ' opened on unit=',luout call openbf(luout,'OUT',luin_tab) ! ! Reset data selection criteria to that for AMSU call set_cloud (idatetime,nfields2d,nlevs,lrc,rcfile,'AMSUA', & f_names,box_size) ! ! Set some variables to start extra loop for AMSU on AQUA nobs2=0 n_mesg_next=1 ierrors2(:)=0 ierrors3(:)=0 n_channels_crtm=15 ! AMSU allocate (crtmBT_M(n_channels_crtm,n_processors)) allocate (lerror2_M(n_processors)) allocate (obs_info_M(obs_ndim4,n_processors)) allocate (obs_values_M(obs_ndim1,obs_ndim2,n_processors)) nchan=n_channels_crtm print *,' ' print ('(3a,i10)'),' 3rd loop Processing subset ',subset1, & ' for date ',idate1 ! ! If the first message in the original file had only the date/time stamp with ! no reports, then in the file to be created, do likewise so that GSI reads ! the file properly. call openmg(luout,subset1,idate1) ! add subset and idate to new file call minimg(luout,0) ! add minutes=00 to idate1 write if (lmsg1) call closmg(luout) ! ! Begin loop over reports do i_obs=1,ierrors(7),n_processors j_obs_last=min (i_obs + n_processors - 1, ierrors(7)) j_obs_span=j_obs_last - i_obs + 1 ! number of obs current set ! ! Read in next set of obs and create profiles do j_obs=i_obs, j_obs_last ! ! Read next report. If no further reports exist in this message, then read ! through messages until next report is found. Reports should always be ! found here since they were written in the previous section of the code, ! so that the abnormal stop should never occur. i_return=ireadsb(luin) if (i_return /= 0) then if (j_obs /=1 ) call closmg(luout) do i=n_mesg_next,ierrors(7) ! look through remaining messages if (ireadmg(luin,subset,idate) == 0) then n_mesg(2)=n_mesg(2)+1 ! count messages in temp file n_mesg_next=n_mesg_next+1 if (ireadsb(luin) == 0) then i_return=1 ! new report exists else i_return=3 ! no reports in message endif else i_return=2 ! no more messages endif if (i_return < 3 ) exit enddo ! if (i_return ==2) then print *,' ' print *,' No more messages readable in third loop' print *,' j_obs, subset, idate =',j_obs,subset,idate stop endif endif ! test on i_return from ireadsb ! ! Read temporary output file that contains simulated AIRS observations ! but extraneous copies of a single real AMSU onservation set call read_write_obs_airs (obs_ndim1,obs_ndim2,obs_ndim3, & obs_ndim4,obs_n_channels,obs_n_data,nobs2, & ierrors2(9),luin,luout,ltest,lsim,bmiss, & obs_info,obs_channels,obs_values,.true.) call fill_prof_info (idatetime,n_errors,ierrors2, & obs_info(n_spot+1:2*n_spot),d_type,prof_info,lerror1) subtype=nint(prof_info(1)) ! ! Interpolate required 2d and 3d fields to obs location if (.not.lerror1) then call interp_fields_for_rad (nlevs,nfields3d,nfields2d, & prof_ndim4,n_spot,nobs2,n_errors,ierrors2, & prof_info,prof_2d,prof_3d,prof_plevs, & .true.,lerror2,ltest) else ! do not process; just flag as error lerror2=.true. endif ! test on .not.lerror1 ! if (lerror2) then ! For invalid profiles, simply copy a previously saved valid profile used for ! AIRS so that the CRTM executes without error, but later set result to ! missing values, as flagged by the lerror2_M array. prof_2d(:)=prof_2d_S(:) prof_3d(:,:)=prof_3d_S(:,:) prof_plevs(:)=prof_plevs_S(:) prof_info(1)=570.0_rkind1 ! indicates AMSUA on AQUA prof_info(2)=obs_info(2) ! copy lat from valid AIRS header prof_info(3)=obs_info(3) ! copy lon from valid AIRS header prof_info(4:n_spot)=prof_info_S(4:n_spot) prof_info(n_spot+1)=-one_k1 ! flag indicates fake profile lcldtop=nlevs prof_info(n_spot+2)=real(lcldtop) pcldtop=prof_2d(1) ! assumes 1st value is ps prof_info(n_spot+3)=pcldtop prof_info(n_spot+3)=zero_k1 ! hours from synoptic time soza=0. else ! For valid profiles, determine observation quality based on cloud or precip ! conditions and count conditions of various types call cloud (nlevs,nfields2d,prof_2d,prof_plevs, & lcldtop,pcldtop) prof_info(n_spot+1)=-2. ! flag indicates correct profile prof_info(n_spot+2)=real(lcldtop) prof_info(n_spot+3)=pcldtop if (lcldtop == nlevs) then cloudy(1,3)=cloudy(1,3)+one_k1 ! cloud free observation else do i=2,n_clouds if ((pcldtop/prof_2d(1) < cloudy(i-1,1)) & .and. (pcldtop/prof_2d(1) >= cloudy(i,1))) then cloudy(i,3)=cloudy(i,3)+one_k1 endif enddo endif soza=obs_info(3*n_spot+6) endif ! ! Copy profile of each obs in set to the array of multiple profiles j_obs_M=j_obs-i_obs+1 prof_info_M(:,j_obs_M)=prof_info(:) prof_plevs_M(:,j_obs_M)=prof_plevs(:) prof_2d_M(:,j_obs_M)=prof_2d(:) prof_3d_M(:,:,j_obs_M)=prof_3d(:,:) lcldtop_M(j_obs_M)=lcldtop pcldtop_M(j_obs_M)=pcldtop soza_M(j_obs_M)=soza lerror2_M(j_obs_M)=lerror2 obs_info_M(:,j_obs_M)=obs_info(:) obs_values_M(:,:,j_obs_M)=obs_values(:,:) ! enddo ! end first loop over j_obs for AMSUA/AQUA ! ! If this is first occurance of subtype, initialize some values if (i_obs == 1) then subtype=nint(prof_info(1)) ! instrument ID # lfirst=.true. ! ! Initialize multiple copies of the CRTM !MMMM (4) This is exactly like a MMMM (2) do j_obs_M=1,1 !MMMM in MPI version, must be =1,j_obs_span error_status = crtm_destroy(channelinfo) if (error_status /= success) then print *,'SIM_OBS_RAD: ***ERROR*** crtm_destroy ', & 'error_status=',error_status endif call InitGetbt_crtm (prof_ndim4,nlevs,nfields2d,nfields3d, & prof_info_M(:,j_obs_M),prof_plevs_M(:,j_obs_M), & prof_2d_M(:,j_obs_M),prof_3d_M(:,:,j_obs_M), & lcldtop_M(j_obs_M),pcldtop_M(j_obs_M), & lfirst,d_type,subtype,crtmBT_M(:,j_obs_M),nchan, & channelinfo,soza_M(j_obs_M),sfc_type_M(j_obs_M)) enddo ! print *, 'Old channelinfo destroyed: new subtype= ', subtype endif ! test on i-obs == 1 ! ! Call CRTM to create sets of observations from sets of profiles lfirst=.false. !MMMM (5) This is exactly like a MMMM (3) do j_obs_M=1,j_obs_span call InitGetbt_crtm (prof_ndim4,nlevs,nfields2d,nfields3d, & prof_info_M(1,j_obs_M),prof_plevs_M(1,j_obs_M), & prof_2d_M(1,j_obs_M),prof_3d_M(1,1,j_obs_M), & lcldtop_M(j_obs_M),pcldtop_M(j_obs_M), & lfirst,d_type,subtype,crtmBT_M(1,j_obs_M),nchan, & channelinfo,soza_M(j_obs_M),sfc_type_M(j_obs_M)) enddo ! ! Write report of simulated obs do j_obs=i_obs, j_obs_last ! 2nd loop over j_obs for AMSUA/AQUA ! ! Copy simulated for obs values at one location from set of multiple obs j_obs_M=j_obs - i_obs + 1 sfc_type=sfc_type_M(j_obs_M) sfc_types(sfc_type,2)=sfc_types(sfc_type,2)+1 prof_info(:)=prof_info_M(:,j_obs_M) obs_info(:)=obs_info_M(:,j_obs_M) obs_values(:,:)=obs_values_M(:,:,j_obs_M) ierrors2(7)=ierrors2(7)+1 ! number of thinned obs. ! ! If errors with AMSU observation profile_info were detected, then set ! to missing value; otherwise copy values to correct locations in array. if (lerror2_M(j_obs_M)) then obs_values(282:281+n_channels_crtm,1)=bmiss else obs_values(282:281+n_channels_crtm,1)= & crtmBT_M(1:n_channels_crtm,j_obs_M) endif ! ! Fill HSB portion of report with missing data obs_values(282+n_channels_crtm:obs_n_channels,1)=bmiss ! call openmb(luout,subset1,idate1) call read_write_obs_airs (obs_ndim1,obs_ndim2,obs_ndim3, & obs_ndim4,obs_n_channels,obs_n_data,ierrors2(7), & ierrors2(9),luin,luout,ltest,lsim,bmiss, & obs_info,obs_channels,obs_values,.false.) ! ! Print sample of data if testing code if (ltest_sample .and. & (mod(ierrors2(7),box_nums(1,2)/4) == 1) ) then ix1=n_spot*n_subtypes+1 ix2=ix1+n_spot2-1 ix3=ix2-ix1+1 call print_sample ('AQUA/AMSUA',ierrors2(7),subtype, & ierrors2(7),prof_ndim4,n_spot,n_spot2, & ix3,n_channels_crtm,sfc_type, & prof_info,obs_info(n_spot+1:2*n_spot), & obs_info(ix1:ix2), & obs_values(282:281+n_channels_crtm,1)) endif ! test on ltest_sample, etc. enddo ! second loop over j_obs for AMSUA/AQUA ! enddo ! loop over all observations in sets of n_processors ! ! Deallocate some arrays used in 3rd pass and close files deallocate (prof_info_S, prof_plevs_S, prof_2d_S, prof_3d_S) deallocate (crtmBT_M, lerror2_M, obs_info_M, obs_values_M) call closmg(luout) ! final close out of messages if (luin /= luin_tab) call closbf(luin_tab) call closbf(luin) call closbf(luout) ! endif ! test on d_type== 'AIRS_' ! ! End 3rd pass through data ! ! ! Print out summary table of number of obs reports processed print *,' ' print *,' SUMMARY TABLE;' print ('(i7,2a)'),nobs,' observation reports read for ',d_type print ('(i7,2a)'),n_mesg(1),' observation message-groups read all types' print ('(i7,2a)'),ierrors(8),' number of empty thinning boxes of', & ' all sub-types' frac=real(ierrors(7))/real(ierrors(8)+ierrors(7)) print ('(f7.4,a)'),frac,' fraction of non-empty boxes' print ('(i7,a)'),ierrors(7),' number of observation reports written out' frac=real(ierrors(7))/real(ierrors(6)) print ('(f7.5,a)'),frac,' fraction of reports written out vs. read in' ! ! Print fractions of simulated observations determined for elevated surfaces cloudy(:,2)=cloudy(:,2)/ierrors(7) print *,' Fractions of simulated observation with surface set as:' print ('(f7.4,a)'),cloudy(1,2),' have surface as actual NR surface' do i=2,n_clouds print ('(f7.4,a,f7.3,a,f6.3)'),cloudy(i,2),' have surface set as ', & cloudy(i-1,1),' > sigma >= ',cloudy(i,1) enddo ! ! Print surface types specified for obs (all IR treated as ocean) print *,' Numbers of profiles for each surface type:' print ('(3(i7,a))'), sfc_types(1,1),'=ocean',sfc_types(2,1),'=ice', & sfc_types(3,1),'=land' ! ! 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(9)+ierrors(10) if (ierrors(n_errors) > 0) then print *,' ' print *,' Summary of bad reports or other errors detected:' print ('(i7,a)'),ierrors(1), & ' observation reports found where timetmax' print ('(i7,a)'),ierrors(3), & ' observation reports found where longitude out of range' print ('(i7,a)'),ierrors(4), & ' observation reports found where latitude out of range' print ('(i7,a)'),ierrors(5), & ' observation reports found but sat or instrument ID not used' print ('(i7,a)'),ierrors(10), & ' errors in creating thinned profiles' print ('(i7,a)'),ierrors(9), & ' errors detected in writing buffer records' endif ! ! Print summary table for AMSUA on AQUA if (d_type == 'AIRS_') then print *,' ' print *,' Summary of AMSUA simulated data on AIRS (AQUA) file' print ('(i7,2a)'),nobs2,' thinned observation reports on temp file' print ('(i7,2a)'),n_mesg(2),' observation message-groups on temp file' print ('(i7,a)'),ierrors2(7),' number of AMSU reports written out' ! ! Print fractions of simulated observations determined for elevated surfaces cloudy(:,3)=cloudy(:,3)/ierrors2(7) print *,' Fractions of simulated observation with surface set as:' print ('(f7.4,a)'),cloudy(1,3),' have surface as actual NR surface' do i=2,n_clouds print ('(f7.4,a,f7.3,a,f6.3)'),cloudy(i,3),' have surface set as ', & cloudy(i-1,1),' > sigma >= ',cloudy(i,1) enddo ! ! Print surface types specified for AMSUA on AQUA print *,' Numbers of profiles for each surface type:' print ('(3(i7,a))'), sfc_types(1,2),'=ocean',sfc_types(2,2),'=ice', & sfc_types(3,2),'=land' ! ! Print out error count summary for AMSUA on AQUA ! ierrors(6)=nobs is number of observation reports read in ierrors2(n_errors)=ierrors2(1)+ierrors2(2)+ierrors2(3)+ierrors2(4) if (ierrors2(n_errors) > 0) then print *,' ' print *,' Summary of bad reports or other errors detected', & ' for AMSUA on AQUA:' print ('(i7,a)'),ierrors2(1), & ' observation reports found where timetmax' print ('(i7,a)'),ierrors2(3), & ' observation reports found where longitude out of range' print ('(i7,a)'),ierrors2(4), & ' observation reports found where latitude out of range' endif ! test on ierrors2(n_errors) ! endif ! test on d_type ! ! Deallocate all arrays deallocate (f_names) deallocate (box_skip, box_nums) deallocate (prof_info, prof_plevs, prof_2d, prof_3d) deallocate (prof_info_M, prof_plevs_M, prof_2d_M, prof_3d_M) deallocate (lcldtop_M, pcldtop_M, soza_M, sfc_type_M) call clean_m_interp ! deallocate grid and field arrays ! print *,' ' print *,'Program completed' ! end program sim_obs_rad ! ! ! X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X ! subroutine print_sample (d_type,num,subtype,ierrors7,prof_ndim4, & n_spot,n_spot2,ndim2,nchan,sfc,prof_info, & obs_info,obs_spot2_info,obs_values) ! ! Print sample of output for comparing runs during testing. ! use m_kinds implicit none character(len=*) :: d_type integer :: num integer :: subtype integer :: ierrors7 integer :: sfc integer :: prof_ndim4, n_spot, n_spot2, ndim2, nchan real(rkind1) :: prof_info(prof_ndim4) real(rkind3) :: obs_info(n_spot) real(rkind3) :: obs_spot2_info(ndim2) real(rkind3) :: obs_values(nchan) ! integer :: n ! print *,' ' print ('(3a,i5,a,i4,a,i6)'),' Sample output: ',d_type,' num_type=', & num,' subtype=',subtype,' num_obs=',ierrors7 print ('(a,i2,a)'),' sfc_type_index=',sfc,' prof_info=' print ('(1p10e13.5)'),prof_info print *,' instrument_obs_info=' print ('(1p10e13.5)'),obs_info if (n_spot2 > 0) then print *,' sat_spot_info=' print ('(1p10e12.5)'),obs_spot2_info endif print *,' brightness temperatures=' do n=1,nchan print ('(i4,1p1e16.7)'),n,obs_values(n) enddo ! end subroutine print_sample