module readsatobs !$$$ module documentation block ! ! module: readsatobs read data from satellite radiance ! diag* files. ! ! prgmmr: whitaker org: esrl/psd date: 2009-02-23 ! ! abstract: read data from satellite radiance diag* files written out ! by GSI forward operator code. ! ! Public Subroutines: ! get_num_satobs: determine the number of observations to read. ! get_satobs_data: read the data and calculate H(x) for ensemble members. ! write_satobs_data: output diag file with spread ! ! Public Variables: ! ! Modules Used: read_diag ! ! program history log: ! 2009-02-23 Initial version. ! 2016-06-03 Collard - Added changes to allow for historical naming conventions ! 2016-11-29 shlyaeva - updated read routine to calculate linearized H(x) ! added write_ozvobs_data to output ensemble spread ! 2017-12-13 shlyaeva - added netcdf diag read/write capability ! ! attributes: ! language: f95 ! !$$$ use kinds, only: r_kind,i_kind,r_single,r_double use read_diag, only: diag_data_fix_list,diag_header_fix_list,diag_header_chan_list, & diag_data_chan_list,diag_data_extra_list,read_radiag_data,read_radiag_header, & diag_data_name_list, open_radiag, close_radiag use params, only: nsats_rad, dsis, sattypes_rad, npefiles, netcdf_diag, & lupd_satbiasc, use_correlated_oberrs, modelspace_vloc implicit none private public :: get_satobs_data, get_num_satobs, write_satobs_data contains ! get number of radiance observations subroutine get_num_satobs(obspath,datestring,num_obs_tot,num_obs_totdiag,id) implicit none character(len=500), intent(in) :: obspath character(len=10), intent(in) :: id, datestring integer(i_kind), intent(out) :: num_obs_tot, num_obs_totdiag if (netcdf_diag) then call get_num_satobs_nc(obspath,datestring,num_obs_tot,num_obs_totdiag,id) else call get_num_satobs_bin(obspath,datestring,num_obs_tot,num_obs_totdiag,id) endif end subroutine get_num_satobs ! get number of radiance observations from binary file subroutine get_num_satobs_bin(obspath,datestring,num_obs_tot,num_obs_totdiag,id) use radinfo, only: iuse_rad,nusis,jpch_rad,npred implicit none character(len=500), intent(in) :: obspath character(len=10), intent(in) :: id, datestring integer(i_kind), intent(out) :: num_obs_tot, num_obs_totdiag character(len=500) obsfile character(len=20) :: sat_type character(len=4) :: pe_name integer(i_kind) iunit, iflag, nsat, n, nkeep, i, jpchstart,indxsat,ipe integer(i_kind) npred_radiag, istatus logical fexist,lretrieval,lverbose,init_pass real(r_kind) :: errorlimit,errorlimit2 type(diag_header_fix_list ) :: header_fix0 type(diag_header_chan_list),allocatable :: header_chan0(:) type(diag_data_fix_list ) :: data_fix0 type(diag_data_chan_list ),allocatable :: data_chan0(:) type(diag_data_extra_list) ,allocatable :: data_extra0(:,:) type(diag_data_name_list) :: data_name0 ! make consistent with screenobs errorlimit=1._r_kind/sqrt(1.e9_r_kind) errorlimit2=1._r_kind/sqrt(1.e-6_r_kind) iunit = 7 lretrieval=.false. print *,'npred = ',npred npred_radiag=npred lverbose=.false. num_obs_tot = 0 num_obs_totdiag = 0 do nsat=1,nsats_rad jpchstart=0 do i=1,jpch_rad write(sat_type,'(a20)') adjustl(dsis(nsat)) ! The following is to sort out some historical naming conventions select case (sat_type(1:4)) case ('airs') sat_type='airs_aqua' case ('iasi') if (index(sat_type,'metop-a') /= 0) sat_type='iasi_metop-a' if (index(sat_type,'metop-b') /= 0) sat_type='iasi_metop-b' if (index(sat_type,'metop-c') /= 0) sat_type='iasi_metop-c' end select if(sat_type == trim(nusis(i)) .and. iuse_rad(i) > 0) then jpchstart=i exit end if end do if(jpchstart == 0) cycle init_pass = .true. peloop: do ipe=0,npefiles write(pe_name,'(i4.4)') ipe if (npefiles .eq. 0) then ! read diag file (concatenated pe* files) obsfile = trim(adjustl(obspath))//"diag_"//trim(sattypes_rad(nsat))//"_ges."//datestring//'_'//trim(adjustl(id)) inquire(file=obsfile,exist=fexist) if (.not. fexist .or. datestring .eq. '0000000000') & obsfile = trim(adjustl(obspath))//"diag_"//trim(sattypes_rad(nsat))//"_ges."//trim(adjustl(id)) else ! read raw, unconcatenated pe* files. obsfile =& trim(adjustl(obspath))//'gsitmp_'//trim(adjustl(id))//'/pe'//pe_name//'.'//trim(sattypes_rad(nsat))//'_01' endif inquire(file=obsfile,exist=fexist) if (.not.fexist) cycle peloop nkeep = 0 call open_radiag(obsfile,iunit,istatus) if (init_pass) then call read_radiag_header(iunit,npred_radiag,lretrieval,header_fix0,header_chan0,data_name0,iflag,lverbose) if (iflag /= 0) exit init_pass = .false. endif do call read_radiag_data(iunit,header_fix0,lretrieval,data_fix0,data_chan0,data_extra0,iflag) if( iflag /= 0 )exit chan: do n=1,header_fix0%nchan num_obs_totdiag = num_obs_totdiag + 1 if(header_chan0(n)%iuse<1) cycle chan indxsat=header_chan0(n)%iochan if(data_chan0(n)%qcmark < 0. .or. data_chan0(n)%errinv < errorlimit & .or. data_chan0(n)%errinv > errorlimit2 & .or. indxsat == 0) cycle chan if (header_fix0%iextra > 0) then if(.not. modelspace_vloc .and. (data_extra0(1,n)%extra <= 0.001_r_kind .or. & data_extra0(1,n)%extra > 1200._r_kind)) cycle chan if(abs(data_chan0(n)%tbobs) > 1.e9_r_kind) cycle chan else if(abs(data_chan0(n)%tbobs) > 1.e9_r_kind) cycle chan endif nkeep = nkeep + 1 end do chan enddo num_obs_tot = num_obs_tot + nkeep call close_radiag(obsfile,iunit) if (ipe .eq. npefiles) then write(6,100) nsat,trim(sattypes_rad(nsat)),num_obs_tot 100 format(2x,i3,2x,a20,2x,'num_obs_tot= ',i9) endif enddo peloop ! ipe enddo ! satellite end subroutine get_num_satobs_bin ! get number of radiance observations from netcdf file subroutine get_num_satobs_nc(obspath,datestring,num_obs_tot,num_obs_totdiag,id) use radinfo, only: iuse_rad,nusis,jpch_rad use nc_diag_read_mod, only: nc_diag_read_get_var use nc_diag_read_mod, only: nc_diag_read_get_dim use nc_diag_read_mod, only: nc_diag_read_init, nc_diag_read_close implicit none character(len=500), intent(in) :: obspath character(len=10), intent(in) :: id, datestring integer(i_kind), intent(out) :: num_obs_tot, num_obs_totdiag character(len=500) obsfile character(len=20) :: sat_type character(len=4) :: pe_name integer(i_kind) iunit, nsat, nobs, nchans, ipe, nkeep, i, jpchstart logical fexist real(r_kind) :: errorlimit,errorlimit2 integer(i_kind), dimension(:), allocatable :: Satinfo_Chan, Use_Flag, chind real(r_single), dimension(:), allocatable :: Pressure, QC_Flag, Inv_Error, Observation ! make consistent with screenobs errorlimit=1._r_kind/sqrt(1.e9_r_kind) errorlimit2=1._r_kind/sqrt(1.e-6_r_kind) iunit = 7 num_obs_tot = 0 num_obs_totdiag = 0 do nsat=1,nsats_rad jpchstart=0 do i=1,jpch_rad write(sat_type,'(a20)') adjustl(dsis(nsat)) ! The following is to sort out some historical naming conventions select case (sat_type(1:4)) case ('airs') sat_type='airs_aqua' case ('iasi') if (index(sat_type,'metop-a') /= 0) sat_type='iasi_metop-a' if (index(sat_type,'metop-b') /= 0) sat_type='iasi_metop-b' if (index(sat_type,'metop-c') /= 0) sat_type='iasi_metop-c' end select if(sat_type == trim(nusis(i)) .and. iuse_rad(i) > 0) then jpchstart=i exit end if end do if(jpchstart == 0) cycle peloop: do ipe=0,npefiles write(pe_name,'(i4.4)') ipe if (npefiles .eq. 0) then ! read diag file (concatenated pe* files) obsfile = trim(adjustl(obspath))//"diag_"//trim(sattypes_rad(nsat))//"_ges."//datestring//'_'//trim(adjustl(id))//".nc4" inquire(file=obsfile,exist=fexist) if (.not. fexist .or. datestring .eq. '0000000000') & obsfile = trim(adjustl(obspath))//"diag_"//trim(sattypes_rad(nsat))//"_ges."//trim(adjustl(id))//".nc4" else ! read raw, unconcatenated pe* files. obsfile =& trim(adjustl(obspath))//'gsitmp_'//trim(adjustl(id))//'/pe'//pe_name//'.'//trim(sattypes_rad(nsat))//'_01'//".nc4" endif inquire(file=obsfile,exist=fexist) if (.not.fexist) cycle peloop nkeep = 0 call nc_diag_read_init(obsfile, iunit) nobs = nc_diag_read_get_dim(iunit,'nobs') if (nobs <= 0) then call nc_diag_read_close(obsfile) cycle peloop endif nchans = nc_diag_read_get_dim(iunit,'nchans') allocate(Satinfo_Chan(nchans), Use_Flag(nchans), Pressure(nobs), QC_Flag(nobs), & Inv_Error(nobs), Observation(nobs), chind(nobs)) call nc_diag_read_get_var(iunit, 'satinfo_chan', Satinfo_Chan) call nc_diag_read_get_var(iunit, 'use_flag', Use_Flag) call nc_diag_read_get_var(iunit, 'Channel_Index', chind) call nc_diag_read_get_var(iunit, 'Press_Max_Weight_Function', Pressure) call nc_diag_read_get_var(iunit, 'QC_Flag', QC_Flag) call nc_diag_read_get_var(iunit, 'Inverse_Observation_Error', Inv_Error) call nc_diag_read_get_var(iunit, 'Observation', Observation) call nc_diag_read_close(obsfile) do i = 1, nobs num_obs_totdiag = num_obs_totdiag + 1 if(Use_Flag(chind(i)) < 1 ) cycle if(QC_Flag(i) < 0. .or. Inv_Error(i) < errorlimit & .or. Inv_Error(i) > errorlimit2 & .or. Satinfo_Chan(chind(i)) == 0) cycle if(.not. modelspace_vloc .and. (Pressure(i) <= 0.001_r_kind .or. & Pressure(i) > 1200._r_kind)) cycle if(abs(Observation(i)) > 1.e9_r_kind) cycle nkeep = nkeep + 1 enddo num_obs_tot = num_obs_tot + nkeep if (ipe .eq. npefiles) then write(6,100) nsat,trim(sattypes_rad(nsat)),num_obs_tot 100 format(2x,i3,2x,a20,2x,'num_obs_tot= ',i9) endif deallocate(Satinfo_Chan, Use_Flag, Pressure, QC_Flag, Inv_Error, Observation, chind) enddo peloop ! ipe enddo ! satellite end subroutine get_num_satobs_nc ! read radiance data subroutine get_satobs_data(obspath, datestring, nobs_max, nobs_maxdiag, hx_mean, hx_mean_nobc, hx, hx_modens, x_obs, x_err, & x_lon, x_lat, x_press, x_time, x_channum, x_errorig, x_type, x_biaspred, x_indx, x_used, id, nanal, nmem) use radinfo, only: npred use params, only: neigv implicit none character*500, intent(in) :: obspath character(len=10), intent(in) :: datestring integer(i_kind), intent(in) :: nobs_max, nobs_maxdiag real(r_single), dimension(nobs_max), intent(out) :: hx_mean,hx_mean_nobc, hx ! hx_modens holds modulated ensemble in ob space (zero size and not referenced if neigv=0) real(r_single), dimension(neigv, nobs_max), intent(out) :: hx_modens real(r_single), dimension(nobs_max), intent(out) :: x_obs real(r_single), dimension(nobs_max), intent(out) :: x_err, x_errorig real(r_single), dimension(nobs_max), intent(out) :: x_lon, x_lat real(r_single), dimension(nobs_max), intent(out) :: x_press, x_time integer(i_kind), dimension(nobs_max), intent(out) :: x_channum, x_indx character(len=20), dimension(nobs_max), intent(out) :: x_type real(r_single), dimension(npred+1,nobs_max), intent(out) :: x_biaspred integer(i_kind), dimension(nobs_maxdiag), intent(out) :: x_used character(len=10), intent(in) :: id integer(i_kind), intent(in) :: nanal,nmem if (netcdf_diag) then call get_satobs_data_nc(obspath, datestring, nobs_max, nobs_maxdiag, hx_mean, hx_mean_nobc, hx, hx_modens, x_obs, x_err, & x_lon, x_lat, x_press, x_time, x_channum, x_errorig, x_type, x_biaspred, x_indx, x_used, id, nanal, nmem) else call get_satobs_data_bin(obspath, datestring, nobs_max, nobs_maxdiag, hx_mean, hx_mean_nobc, hx, hx_modens, x_obs, x_err, & x_lon, x_lat, x_press, x_time, x_channum, x_errorig, x_type, x_biaspred, x_indx, x_used, id, nanal, nmem) endif end subroutine get_satobs_data ! read radiance data from binary file subroutine get_satobs_data_bin(obspath, datestring, nobs_max, nobs_maxdiag, hx_mean, hx_mean_nobc, hx, hx_modens, x_obs, x_err, & x_lon, x_lat, x_press, x_time, x_channum, x_errorig, x_type, x_biaspred, x_indx, x_used, id, nanal, nmem) use radinfo, only: iuse_rad,nusis,jpch_rad,npred,adp_anglebc,emiss_bc use params, only: nanals, lobsdiag_forenkf, neigv, vlocal_evecs use statevec, only: state_d use constants, only: deg2rad, zero use mpisetup, only: nproc, mpi_wtime use observer_enkf, only: calc_linhx, calc_linhx_modens, setup_linhx use sparsearr, only: sparr, init_raggedarr, raggedarr implicit none character*500, intent(in) :: obspath character(len=10), intent(in) :: datestring type(raggedarr) :: hxpert integer(i_kind), intent(in) :: nobs_max, nobs_maxdiag real(r_single), dimension(nobs_max), intent(out) :: hx_mean,hx_mean_nobc, hx ! hx_modens holds modulated ensemble in ob space (zero size and not referenced if neigv=0) real(r_single), dimension(neigv, nobs_max), intent(out) :: hx_modens real(r_single), dimension(nobs_max), intent(out) :: x_obs real(r_single), dimension(nobs_max), intent(out) :: x_err, x_errorig real(r_single), dimension(nobs_max), intent(out) :: x_lon, x_lat real(r_single), dimension(nobs_max), intent(out) :: x_press, x_time integer(i_kind), dimension(nobs_max), intent(out) :: x_channum, x_indx character(len=20), dimension(nobs_max), intent(out) :: x_type real(r_single), dimension(npred+1,nobs_max), intent(out) :: x_biaspred integer(i_kind), dimension(nobs_maxdiag), intent(out) :: x_used character(len=10), intent(in) :: id integer(i_kind), intent(in) :: nanal, nmem character*500 obsfile, obsfile2 character(len=10) :: id2 character(len=4) pe_name character(len=20) :: sat_type integer(i_kind) iunit, iflag,nob,nobdiag, n,nsat,ipe,i,jpchstart,indxsat,nn integer(i_kind) iunit2, iflag2, istatus integer(i_kind) npred_radiag logical fexist,lretrieval,lverbose,init_pass logical twofiles,fexist2,init_pass2 real(r_kind) :: errorlimit,errorlimit2 real(r_double) t1,t2,tsum,tsum2 integer(i_kind) :: ix, iy, it, ixp, iyp, itp real(r_kind) :: delx, dely, delxp, delyp, delt, deltp real(r_single) :: rlat,rlon,rtim,rlat_prev,rlon_prev,rtim_prev,eps type(diag_header_fix_list ) :: header_fix, header_fix2 type(diag_header_chan_list),allocatable :: header_chan(:), header_chan2(:) type(diag_data_fix_list ) :: data_fix, data_fix2 type(diag_data_chan_list ),allocatable :: data_chan(:), data_chan2(:) type(diag_data_extra_list) ,allocatable :: data_extra(:,:), data_extra2(:,:) type(diag_data_name_list) :: data_name, data_name2 ! make consistent with screenobs errorlimit=1._r_kind/sqrt(1.e9_r_kind) errorlimit2=1._r_kind/sqrt(1.e-6_r_kind) eps = 1.e-3 tsum = 0; tsum2 = 0 iunit = 7 iunit2 = 17 lretrieval=.false. npred_radiag=npred lverbose=.false. twofiles = (.not. lobsdiag_forenkf) .and. (nanal <= nanals) id2 = 'ensmean' if (nanal <= nanals) then write(id2,'(a3,(i3.3))') 'mem',nanal endif hx = zero nob = 0 rlat_prev = -1.e30; rlon_prev=-1.e30; rtim_prev = -1.e30 nobdiag = 0 x_used = 0 do nsat=1,nsats_rad jpchstart=0 do i=1,jpch_rad write(sat_type,'(a20)') adjustl(dsis(nsat)) ! The following is to sort out some historical naming conventions select case (sat_type(1:4)) case ('airs') sat_type='airs_aqua' case ('iasi') if (index(sat_type,'metop-a') /= 0) sat_type='iasi_metop-a' if (index(sat_type,'metop-b') /= 0) sat_type='iasi_metop-b' if (index(sat_type,'metop-c') /= 0) sat_type='iasi_metop-c' end select if(sat_type == trim(nusis(i)) .and. iuse_rad(i) > 0) then jpchstart = i exit end if end do if(jpchstart == 0) cycle init_pass = .true.; init_pass2 = .true. peloop: do ipe=0,npefiles write(pe_name,'(i4.4)') ipe if (npefiles .eq. 0) then ! read diag file (concatenated pe* files) obsfile = trim(adjustl(obspath))//"diag_"//trim(sattypes_rad(nsat))//"_ges."//datestring//'_'//trim(adjustl(id)) inquire(file=obsfile,exist=fexist) if (.not. fexist .or. datestring .eq. '0000000000') & obsfile = trim(adjustl(obspath))//"diag_"//trim(sattypes_rad(nsat))//"_ges."//trim(adjustl(id)) else ! read raw, unconcatenated pe* files. obsfile =& trim(adjustl(obspath))//'gsitmp_'//trim(adjustl(id))//'/pe'//pe_name//'.'//trim(sattypes_rad(nsat))//'_01' endif inquire(file=obsfile,exist=fexist) if(.not.fexist) cycle peloop call open_radiag(obsfile,iunit,istatus) if (init_pass) then call read_radiag_header(iunit,npred_radiag,lretrieval,header_fix,header_chan,data_name,iflag,lverbose) if( iflag /= 0 ) exit init_pass = .false. endif if(twofiles)then if (npefiles .eq. 0) then ! read diag file (concatenated pe* files) obsfile2 = trim(adjustl(obspath))//"diag_"//trim(sattypes_rad(nsat))//"_ges."//datestring//'_'//trim(adjustl(id2)) inquire(file=obsfile2,exist=fexist2) if (.not. fexist2 .or. datestring .eq. '0000000000') & obsfile2 = trim(adjustl(obspath))//"diag_"//trim(sattypes_rad(nsat))//"_ges."//trim(adjustl(id2)) else ! read raw, unconcatenated pe* files. obsfile2 =& trim(adjustl(obspath))//'gsitmp_'//trim(adjustl(id2))//'/pe'//pe_name//'.'//trim(sattypes_rad(nsat))//'_01' endif call open_radiag(obsfile2, iunit2,istatus) if (init_pass2) then call read_radiag_header(iunit2,npred_radiag,lretrieval,header_fix2,header_chan2,data_name2,iflag2,lverbose) init_pass2 = .false. endif end if do t1 = mpi_wtime() call read_radiag_data(iunit,header_fix,lretrieval,data_fix,data_chan,data_extra,iflag ) t2 = mpi_wtime() tsum2 = tsum2 + t2-t1 if( iflag /= 0 ) then exit end if if(twofiles)then call read_radiag_data(iunit2,header_fix2,lretrieval,data_fix2,data_chan2,data_extra2,iflag2 ) if( header_fix%nchan /= header_fix2%nchan .or. abs(data_fix%lat-data_fix2%lat) .gt. 1.e-5 .or. & abs(data_fix%lon-data_fix2%lon) .gt. 1.e-5 .or. abs(data_fix%obstime-data_fix2%obstime) .gt. 1.e-5) then write(6,*) 'inconsistent files',trim(obsfile2) write(6,*) 'nchan',header_fix%nchan,header_fix2%nchan write(6,*) 'lat',data_fix%lat,data_fix2%lat write(6,*) 'lon',data_fix%lon,data_fix2%lon write(6,*) 'obstim',data_fix%obstime,data_fix2%obstime call stop2(-99) end if end if chan:do n=1,header_fix%nchan nobdiag = nobdiag + 1 if(header_chan(n)%iuse<1) cycle chan indxsat=header_chan(n)%iochan if(data_chan(n)%qcmark < 0. .or. data_chan(n)%errinv < errorlimit & .or. data_chan(n)%errinv > errorlimit2 & .or. indxsat == 0) cycle chan if (header_fix%iextra > 0) then if(.not. modelspace_vloc .and. (data_extra(1,n)%extra <= 0.001_r_kind .or. & data_extra(1,n)%extra > 1200._r_kind)) cycle chan if(abs(data_chan(n)%tbobs) > 1.e9_r_kind) cycle chan else if(abs(data_chan(n)%tbobs) > 1.e9_r_kind) cycle chan endif nob = nob + 1 x_used(nobdiag) = 1 if (nob > nobs_max) then print *,'warning: exceeding array bounds in readinfo_from_file',& nob,nobs_max end if x_type(nob)= sat_type x_channum(nob) = n x_indx(nob) = indxsat x_lon(nob) = data_fix%lon x_lat(nob) = data_fix%lat x_time(nob) = data_fix%obstime x_obs(nob) = data_chan(n)%tbobs ! bias corrected Hx hx_mean(nob) = x_obs(nob) - data_chan(n)%omgbc ! un-bias corrected Hx hx_mean_nobc(nob) = x_obs(nob) - data_chan(n)%omgnbc if (nanal <= nanals) then ! read full Hx if (.not. lobsdiag_forenkf) then hx(nob) = x_obs(nob) - data_chan2(n)%omgbc ! run linearized Hx else rlat = x_lat(nob)*deg2rad rlon = x_lon(nob)*deg2rad rtim = x_time(nob) if (nob > 1) then rlat_prev = x_lat(nob-1)*deg2rad rlon_prev = x_lon(nob-1)*deg2rad rtim_prev = x_time(nob-1) endif if (abs(rlat-rlat_prev) > eps .or. & abs(rlon-rlon_prev) > eps .or. & abs(rtim-rtim_prev) > eps) then call setup_linhx(rlat,rlon,rtim, & ix, delx, ixp, delxp, iy, dely, & iyp, delyp, it, delt, itp, deltp) endif call init_raggedarr(hxpert, data_chan(n)%dhx_dx%nnz) call calc_linhx(hx_mean(nob), state_d(:,:,:,nmem), & data_chan(n)%dhx_dx, hxpert, hx(nob), & ix, delx, ixp, delxp, iy, dely, & iyp, delyp, it, delt, itp, deltp) ! compute modulated ensemble in obs space if (neigv>0) call calc_linhx_modens(hx_mean(nob),data_chan(n)%dhx_dx,hxpert,hx_modens(:,nob),vlocal_evecs) t2 = mpi_wtime() tsum = tsum + t2-t1 endif endif ! data_chan%errinv is inverse error variance. x_errorig(nob) = header_chan(n)%varch**2 x_err(nob) = (1._r_kind/data_chan(n)%errinv)**2 if (header_fix%iextra > 0) then x_press(nob) = data_extra(1,n)%extra else x_press(nob) = 99999 endif !! DTK: **NOTE** !! The bifix term will need to be expanded if/when the GSI/GDAS goes to using !! a higher polynomial version of the angle dependent bias correction (if !! and when it is moved into part of the varbc) ! from radinfo: radiance bias correction terms are as follows: ! pred(1,:) = global offset ! pred(2,:) = zenith angle predictor, is not used and set to zero now ! pred(3,:) = cloud liquid water predictor for clear-sky microwave radiance assimilation ! pred(4,:) = square of temperature laps rate predictor ! pred(5,:) = temperature laps rate predictor ! pred(6,:) = cosinusoidal predictor for SSMI/S ascending/descending bias ! pred(7,:) = sinusoidal predictor for SSMI/S ! pred(8,:) = emissivity sensitivity predictor for land/sea differences ! pred(9,:) = fourth order polynomial of angle bias correction ! pred(10,:) = third order polynomial of angle bias correction ! pred(11,:) = second order polynomial of angle bias correction ! pred(12,:) = first order polynomial of angle bias correction if (lupd_satbiasc) then ! bias predictors only used if lupd_satbiasc=T x_biaspred(1,nob) = data_chan(n)%bicons ! constant bias correction x_biaspred(2,nob) = data_chan(n)%biang ! scan angle bias correction x_biaspred(3,nob) = data_chan(n)%biclw ! CLW bias correction x_biaspred(4,nob) = data_chan(n)%bilap2 ! square lapse rate bias corr x_biaspred(5,nob) = data_chan(n)%bilap ! lapse rate bias correction x_biaspred(6,nob) = data_chan(n)%bicos ! node*cos(lat) bias correction for SSMIS x_biaspred(7,nob) = data_chan(n)%bisin ! sin(lat) bias correction for SSMIS if (emiss_bc) then x_biaspred(8,nob) = data_chan(n)%biemis nn = 9 else nn = 8 endif if (adp_anglebc) then x_biaspred(nn ,nob) = data_chan(n)%bifix(1) ! 4th order scan angle (predictor) x_biaspred(nn+1,nob) = data_chan(n)%bifix(2) ! 3rd order scan angle (predictor) x_biaspred(nn+2,nob) = data_chan(n)%bifix(3) ! 2nd order scan angle (predictor) x_biaspred(nn+3,nob) = data_chan(n)%bifix(4) ! 1st order scan angle (predictor) endif x_biaspred(npred+1,nob) = data_chan(n)%bifix(1) ! fixed angle dependent bias else x_biaspred(:,nob) =zero ! lupd_satbiasc=F, don't need bias predictors endif enddo chan enddo call close_radiag(obsfile,iunit) if (twofiles) call close_radiag(obsfile2,iunit2) enddo peloop ! ipe enddo ! satellite if (nanal == nanals .and. lobsdiag_forenkf) print *,'time in calc_linhx for sat obs on proc',nproc,' = ',tsum if (nanal == nanals) print *,'time in read_raddiag_data for sat obs on proc',nproc,' = ',tsum2 if (nob /= nobs_max) then print *,'number of obs not what expected in get_satobs_data',nob,nobs_max call stop2(92) end if if (nobdiag /= nobs_maxdiag) then print *,'number of total diag obs not what expected in get_satobs_data',nobdiag,nobs_maxdiag call stop2(92) end if end subroutine get_satobs_data_bin ! read radiance data from netcdf file subroutine get_satobs_data_nc(obspath, datestring, nobs_max, nobs_maxdiag, hx_mean, hx_mean_nobc, hx, hx_modens, x_obs, x_err, & x_lon, x_lat, x_press, x_time, x_channum, x_errorig, x_type, x_biaspred, x_indx, x_used, id, nanal, nmem) use nc_diag_read_mod, only: nc_diag_read_get_var use nc_diag_read_mod, only: nc_diag_read_get_dim, nc_diag_read_get_global_attr use nc_diag_read_mod, only: nc_diag_read_init, nc_diag_read_close use radinfo, only: iuse_rad,nusis,jpch_rad,npred,adp_anglebc,emiss_bc use params, only: nanals, lobsdiag_forenkf, neigv, vlocal_evecs use statevec, only: state_d use constants, only: deg2rad, zero use mpisetup, only: nproc, mpi_wtime use observer_enkf, only: calc_linhx, calc_linhx_modens, setup_linhx use sparsearr, only: sparr, assignment(=), delete, sparr2, new, & init_raggedarr, raggedarr implicit none character*500, intent(in) :: obspath character(len=10), intent(in) :: datestring integer(i_kind), intent(in) :: nobs_max, nobs_maxdiag real(r_single), dimension(nobs_max), intent(out) :: hx_mean,hx_mean_nobc, hx ! hx_modens holds modulated ensemble in ob space (zero size and not referenced if neigv=0) real(r_single), dimension(neigv,nobs_max), intent(out) :: hx_modens real(r_single), dimension(nobs_max), intent(out) :: x_obs real(r_single), dimension(nobs_max), intent(out) :: x_err, x_errorig real(r_single), dimension(nobs_max), intent(out) :: x_lon, x_lat real(r_single), dimension(nobs_max), intent(out) :: x_press, x_time integer(i_kind), dimension(nobs_max), intent(out) :: x_channum, x_indx character(len=20), dimension(nobs_max), intent(out) :: x_type real(r_single), dimension(npred+1,nobs_max), intent(out) :: x_biaspred integer(i_kind), dimension(nobs_maxdiag), intent(out) :: x_used character(len=10), intent(in) :: id integer(i_kind), intent(in) :: nanal, nmem character*500 obsfile, obsfile2 character(len=10) :: id2 character(len=4) pe_name character(len=20) :: sat_type integer(i_kind) iunit, nobs, nobdiag, i, nsat, ipe, jpchstart, nchans integer(i_kind) iunit2, nob, nobs2, nnz, nind, nn integer(i_kind) npred_radiag, angord logical fexist logical twofiles,fexist2 real(r_kind) :: errorlimit,errorlimit2 real(r_double) t1,t2,tsum,tsum2 real(r_single) :: rlat,rlon,rtim,rlat_prev,rlon_prev,rtim_prev,eps type(sparr2) :: dhx_dx_read type(sparr) :: dhx_dx type(raggedarr) :: hxpert integer(i_kind), dimension(:), allocatable :: Satinfo_Chan, Use_Flag, chind, chaninfoidx real(r_kind), dimension(:), allocatable :: error_variance real(r_single), dimension(:), allocatable :: Pressure, QC_Flag, Inv_Error, Inv_Error_scaled, & Observation, Observation_scaled real(r_single), dimension(:), allocatable :: Latitude, Longitude, Time real(r_single), dimension(:), allocatable :: Obs_Minus_Forecast_adjusted real(r_single), dimension(:), allocatable :: Obs_Minus_Forecast_adjusted_scaled real(r_single), dimension(:), allocatable :: Obs_Minus_Forecast_adjusted_scaled2 real(r_single), dimension(:), allocatable :: Obs_Minus_Forecast_unadjusted, Obs_Minus_Forecast_adjusted2 integer(i_kind), allocatable, dimension (:,:) :: Observation_Operator_Jacobian_stind integer(i_kind), allocatable, dimension (:,:) :: Observation_Operator_Jacobian_endind real(r_single), allocatable, dimension (:,:) :: Observation_Operator_Jacobian_val real(r_single), dimension(:), allocatable :: BC_Fixed_Scan_Position, BCPred_Constant, BCPred_Scan_Angle real(r_single), dimension(:), allocatable :: BCPred_Cloud_Liquid_Water, BCPred_Lapse_Rate_Squared, BCPred_Lapse_Rate real(r_single), dimension(:), allocatable :: BCPred_Cosine_Latitude_times_Node, BCPred_Sine_Latitude real(r_single), dimension(:), allocatable :: BCPred_Emissivity real(r_single), allocatable, dimension (:,:) :: BCPred_angord integer(i_kind) :: ix, iy, it, ixp, iyp, itp, nprof real(r_kind) :: delx, dely, delxp, delyp, delt, deltp ! make consistent with screenobs errorlimit=1._r_kind/sqrt(1.e9_r_kind) errorlimit2=1._r_kind/sqrt(1.e-6_r_kind) eps = 1.e-3 tsum = 0; tsum2 = 0 npred_radiag=npred twofiles = (.not. lobsdiag_forenkf) .and. (nanal <= nanals) id2 = 'ensmean' if (nanal <= nanals) then write(id2,'(a3,(i3.3))') 'mem',nanal endif hx = zero nob = 0 rlat_prev = huge(rlat); rlon_prev=huge(rlon); rtim_prev = huge(rtim) nobdiag = 0 x_used = 0 nprof = 0 do nsat=1,nsats_rad jpchstart=0 do i=1,jpch_rad write(sat_type,'(a20)') adjustl(dsis(nsat)) ! The following is to sort out some historical naming conventions select case (sat_type(1:4)) case ('airs') sat_type='airs_aqua' case ('iasi') if (index(sat_type,'metop-a') /= 0) sat_type='iasi_metop-a' if (index(sat_type,'metop-b') /= 0) sat_type='iasi_metop-b' if (index(sat_type,'metop-c') /= 0) sat_type='iasi_metop-c' end select if(sat_type == trim(nusis(i)) .and. iuse_rad(i) > 0) then jpchstart = i exit end if end do if(jpchstart == 0) cycle peloop: do ipe=0,npefiles write(pe_name,'(i4.4)') ipe if (npefiles .eq. 0) then ! read diag file (concatenated pe* files) obsfile = trim(adjustl(obspath))//"diag_"//trim(sattypes_rad(nsat))//"_ges."//datestring//'_'//trim(adjustl(id))//".nc4" inquire(file=obsfile,exist=fexist) if (.not. fexist .or. datestring .eq. '0000000000') & obsfile = trim(adjustl(obspath))//"diag_"//trim(sattypes_rad(nsat))//"_ges."//trim(adjustl(id))//".nc4" else ! read raw, unconcatenated pe* files. obsfile =& trim(adjustl(obspath))//'gsitmp_'//trim(adjustl(id))//'/pe'//pe_name//'.'//trim(sattypes_rad(nsat))//'_01.nc4' endif inquire(file=obsfile,exist=fexist) if(.not.fexist) cycle peloop t1 = mpi_wtime() call nc_diag_read_init(obsfile, iunit) nobs = nc_diag_read_get_dim(iunit,'nobs') if (nobs <= 0) then call nc_diag_read_close(obsfile) cycle peloop endif nchans = nc_diag_read_get_dim(iunit,'nchans') allocate(Satinfo_Chan(nchans), Use_Flag(nchans), error_variance(nchans)) allocate(Pressure(nobs), QC_Flag(nobs), Inv_Error(nobs), Latitude(nobs), & Longitude(nobs), Time(nobs), Observation(nobs), chind(nobs), & Obs_Minus_Forecast_unadjusted(nobs), Obs_Minus_Forecast_adjusted(nobs)) call nc_diag_read_get_var(iunit, 'satinfo_chan', Satinfo_Chan) call nc_diag_read_get_var(iunit, 'use_flag', Use_Flag) call nc_diag_read_get_var(iunit, 'error_variance', error_variance) call nc_diag_read_get_var(iunit, 'chaninfoidx', chaninfoidx) call nc_diag_read_get_var(iunit, 'Channel_Index', chind) call nc_diag_read_get_var(iunit, 'Press_Max_Weight_Function', Pressure) call nc_diag_read_get_var(iunit, 'QC_Flag', QC_Flag) call nc_diag_read_get_var(iunit, 'Inverse_Observation_Error', Inv_Error) call nc_diag_read_get_var(iunit, 'Latitude', Latitude) call nc_diag_read_get_var(iunit, 'Longitude', Longitude) call nc_diag_read_get_var(iunit, 'Obs_Time', Time) call nc_diag_read_get_var(iunit, 'Observation', Observation) call nc_diag_read_get_var(iunit, 'Obs_Minus_Forecast_unadjusted', Obs_Minus_Forecast_unadjusted) call nc_diag_read_get_var(iunit, 'Obs_Minus_Forecast_adjusted', Obs_Minus_Forecast_adjusted) if (use_correlated_oberrs) then call nc_diag_read_get_var(iunit, 'Observation_scaled', Observation_scaled) call nc_diag_read_get_var(iunit, 'Obs_Minus_Forecast_adjusted_scaled', & Obs_Minus_Forecast_adjusted_scaled) call nc_diag_read_get_var(iunit, 'Inverse_Observation_Error_scaled', & Inv_Error_scaled) endif if (lupd_satbiasc) then ! bias predictors only needed if lupd_satbiasc=T allocate(BC_Fixed_Scan_Position(nobs), BCPred_Constant(nobs), BCPred_Scan_Angle(nobs), & BCPred_Cloud_Liquid_Water(nobs), BCPred_Lapse_Rate_Squared(nobs), & BCPred_Lapse_Rate(nobs)) call nc_diag_read_get_var(iunit, 'BC_Fixed_Scan_Position', BC_Fixed_Scan_Position) call nc_diag_read_get_var(iunit, 'BCPred_Constant', BCPred_Constant) call nc_diag_read_get_var(iunit, 'BCPred_Scan_Angle', BCPred_Scan_Angle) call nc_diag_read_get_var(iunit, 'BCPred_Cloud_Liquid_Water', BCPred_Cloud_Liquid_Water) call nc_diag_read_get_var(iunit, 'BCPred_Lapse_Rate_Squared', BCPred_Lapse_Rate_Squared) call nc_diag_read_get_var(iunit, 'BCPred_Lapse_Rate', BCPred_Lapse_Rate) allocate(BCPred_Cosine_Latitude_times_Node(nobs), BCPred_Sine_Latitude(nobs)) call nc_diag_read_get_var(iunit, 'BCPred_Cosine_Latitude_times_Node', BCPred_Cosine_Latitude_times_Node) call nc_diag_read_get_var(iunit, 'BCPred_Sine_Latitude', BCPred_Sine_Latitude) if (emiss_bc) then allocate(BCPred_Emissivity(nobs)) call nc_diag_read_get_var(iunit, 'BCPred_Emissivity', BCPred_Emissivity) endif if (adp_anglebc) then call nc_diag_read_get_global_attr(iunit, "angord", angord) allocate(BCPred_angord(angord, nobs)) call nc_diag_read_get_var(iunit, 'BCPred_angord', BCPred_angord) endif endif ! lupd_satbiasc=T, read bias predictors if (lobsdiag_forenkf) then call nc_diag_read_get_global_attr(iunit, "jac_nnz", nnz) call nc_diag_read_get_global_attr(iunit, "jac_nind", nind) allocate(Observation_Operator_Jacobian_stind(nind, nobs)) allocate(Observation_Operator_Jacobian_endind(nind, nobs)) allocate(Observation_Operator_Jacobian_val(nnz, nobs)) call nc_diag_read_get_var(iunit, 'Observation_Operator_Jacobian_stind', Observation_Operator_Jacobian_stind) call nc_diag_read_get_var(iunit, 'Observation_Operator_Jacobian_endind', Observation_Operator_Jacobian_endind) call nc_diag_read_get_var(iunit, 'Observation_Operator_Jacobian_val', Observation_Operator_Jacobian_val) endif call nc_diag_read_close(obsfile) t2 = mpi_wtime() tsum2 = tsum2 + t2-t1 if(twofiles)then if (npefiles .eq. 0) then ! read diag file (concatenated pe* files) obsfile2 = trim(adjustl(obspath))//"diag_"//trim(sattypes_rad(nsat))//"_ges."//datestring//'_'//trim(adjustl(id2))//".nc4" inquire(file=obsfile2,exist=fexist2) if (.not. fexist2 .or. datestring .eq. '0000000000') & obsfile2 = trim(adjustl(obspath))//"diag_"//trim(sattypes_rad(nsat))//"_ges."//trim(adjustl(id2))//".nc4" else ! read raw, unconcatenated pe* files. obsfile2 =& trim(adjustl(obspath))//'gsitmp_'//trim(adjustl(id2))//'/pe'//pe_name//'.'//trim(sattypes_rad(nsat))//'_01'//".nc4" endif call nc_diag_read_init(obsfile2, iunit2) nobs2 = nc_diag_read_get_dim(iunit2,'nobs') if (nobs2 /= nobs) print *, nanal, trim(obsfile), nobs, nobs2 if (use_correlated_oberrs) then allocate(Obs_Minus_Forecast_adjusted_scaled2(nobs)) call nc_diag_read_get_var(iunit2, 'Obs_Minus_Forecast_adjusted_scaled', & Obs_Minus_Forecast_adjusted_scaled2) else allocate(Obs_Minus_Forecast_adjusted2(nobs)) call nc_diag_read_get_var(iunit2, 'Obs_Minus_Forecast_adjusted', & Obs_Minus_Forecast_adjusted2) endif call nc_diag_read_close(obsfile2) end if do i = 1, nobs nobdiag = nobdiag + 1 if (Use_Flag(chind(i)) < 1) cycle if (QC_Flag(i) < 0. .or. Inv_Error(i) < errorlimit & .or. Inv_Error(i) > errorlimit2 & .or. Satinfo_Chan(chind(i)) == 0) cycle if (.not. modelspace_vloc .and. (Pressure(i) <= 0.001_r_kind .or. & Pressure(i) > 1200._r_kind)) cycle if (abs(Observation(i)) > 1.e9_r_kind) cycle nob = nob + 1 x_used(nobdiag) = 1 x_type(nob)= sat_type x_channum(nob) = chaninfoidx(chind(i)) x_indx(nob) = Satinfo_Chan(chind(i)) x_lon(nob) = Longitude(i) x_lat(nob) = Latitude(i) x_time(nob) = Time(i) ! bias corrected Hx if (use_correlated_oberrs) then x_obs(nob) = Observation_scaled(i) hx_mean(nob) = x_obs(nob) - Obs_Minus_Forecast_adjusted_scaled(i) else x_obs(nob) = Observation(i) hx_mean(nob) = x_obs(nob) - Obs_Minus_Forecast_adjusted(i) endif ! un-bias corrected Hx hx_mean_nobc(nob) = x_obs(nob) - Obs_Minus_Forecast_unadjusted(i) if (nanal <= nanals) then ! read full Hx ! use_correlated_oberrs implies lobsdiag_forenkf if (.not. lobsdiag_forenkf) then hx(nob) = x_obs(nob) - Obs_Minus_Forecast_adjusted2(i) ! run linearized Hx else call new(dhx_dx_read, nnz, nind) dhx_dx_read%st_ind = Observation_Operator_Jacobian_stind(:,i) dhx_dx_read%end_ind = Observation_Operator_Jacobian_endind(:,i) dhx_dx_read%val = Observation_Operator_Jacobian_val(:,i) dhx_dx = dhx_dx_read call init_raggedarr(hxpert, dhx_dx%nnz) t1 = mpi_wtime() rlat = x_lat(nob)*deg2rad rlon = x_lon(nob)*deg2rad rtim = x_time(nob) if (nob > 1) then rlat_prev = x_lat(nob-1)*deg2rad rlon_prev = x_lon(nob-1)*deg2rad rtim_prev = x_time(nob-1) endif if (abs(rlat-rlat_prev) > eps .or. & abs(rlon-rlon_prev) > eps .or. & abs(rtim-rtim_prev) > eps) then call setup_linhx(rlat,rlon,rtim, & ix, delx, ixp, delxp, iy, dely, & iyp, delyp, it, delt, itp, deltp) else nprof = nprof + 1 endif ! note: bias corrected mean added here, but removed later from hx call calc_linhx(hx_mean(nob), state_d(:,:,:,nmem), & dhx_dx, hxpert, hx(nob), & ix, delx, ixp, delxp, iy, dely, & iyp, delyp, it, delt, itp, deltp) ! compute modulated ensemble in obs space if (neigv>0) call calc_linhx_modens(hx_mean(nob),dhx_dx,hxpert,hx_modens(:,nob),vlocal_evecs) t2 = mpi_wtime() tsum = tsum + t2-t1 call delete(dhx_dx) call delete(dhx_dx_read) endif endif x_errorig(nob) = error_variance(chind(i))**2 if (use_correlated_oberrs) then x_err(nob) = (1._r_kind/Inv_Error_scaled(i))**2 else x_err(nob) = (1._r_kind/Inv_Error(i))**2 endif x_press(nob) = Pressure(i) ! DTK: **NOTE** ! The bifix term will need to be expanded if/when the GSI/GDAS goes to using ! a higher polynomial version of the angle dependent bias correction (if ! and when it is moved into part of the varbc) ! from radinfo: radiance bias correction terms are as follows: ! pred(1,:) = global offset ! pred(2,:) = zenith angle predictor, is not used and set to zero now ! pred(3,:) = cloud liquid water predictor for clear-sky microwave radiance assimilation ! pred(4,:) = square of temperature laps rate predictor ! pred(5,:) = temperature laps rate predictor ! pred(6,:) = cosinusoidal predictor for SSMI/S ascending/descending bias ! pred(7,:) = sinusoidal predictor for SSMI/S ! pred(8,:) = emissivity sensitivity predictor for land/sea differences ! pred(9,:) = fourth order polynomial of angle bias correction ! pred(10,:) = third order polynomial of angle bias correction ! pred(11,:) = second order polynomial of angle bias correction ! pred(12,:) = first order polynomial of angle bias correction if (lupd_satbiasc) then ! bias predictors only used if lupd_satbiasc=T x_biaspred(1,nob) = BCPred_Constant(i) ! global offset x_biaspred(2,nob) = BCPred_Scan_Angle(i) ! zenith angle predictor, not used x_biaspred(3,nob) = BCPred_Cloud_Liquid_Water(i) ! CLW bias correction x_biaspred(4,nob) = BCPred_Lapse_Rate_Squared(i) ! square lapse rate bias corr x_biaspred(5,nob) = BCPred_Lapse_Rate(i) ! lapse rate bias correction x_biaspred(6,nob) = BCPred_Cosine_Latitude_times_Node(i) ! node*cos(lat) bias correction for SSMIS x_biaspred(7,nob) = BCPred_Sine_Latitude(i) ! sin(lat) bias correction for SSMIS if (emiss_bc) then x_biaspred(8,nob) = BCPred_Emissivity(i) nn = 9 else nn = 8 endif if (adp_anglebc) then x_biaspred(nn ,nob) = BCPred_angord(1,i) ! 4th order scan angle (predictor) x_biaspred(nn+1,nob) = BCPred_angord(2,i) ! 3rd order scan angle (predictor) x_biaspred(nn+2,nob) = BCPred_angord(3,i) ! 2nd order scan angle (predictor) x_biaspred(nn+3,nob) = BCPred_angord(4,i) ! 1st order scan angle (predictor) endif ! total angle dependent bias correction (sum of four terms) x_biaspred(npred+1,nob) = BC_Fixed_Scan_Position(i) else x_biaspred(:,nob)=zero ! lupd_biaspredc=F, don't need bias predictors endif enddo deallocate(Satinfo_Chan, Use_Flag, error_variance, chaninfoidx) deallocate(Pressure, QC_Flag, Inv_Error, Latitude, Longitude, Time, & Observation, chind, Obs_Minus_Forecast_unadjusted, & Obs_Minus_Forecast_adjusted) if (use_correlated_oberrs) then deallocate(Obs_Minus_Forecast_adjusted_scaled,Inv_Error_scaled,& Observation_scaled) endif if (lupd_satbiasc) then ! bias predictors only used if lupd_satbiasc=T deallocate(BC_Fixed_Scan_Position, BCPred_Constant, BCPred_Scan_Angle, & BCPred_Cloud_Liquid_Water, BCPred_Lapse_Rate_Squared, & BCPred_Lapse_Rate) deallocate(BCPred_Cosine_Latitude_times_Node, BCPred_Sine_Latitude) if (emiss_bc) deallocate(BCPred_Emissivity) if (adp_anglebc) deallocate(BCPred_angord) endif if (twofiles) then deallocate(Obs_Minus_Forecast_adjusted2) if (use_correlated_oberrs) deallocate(Obs_Minus_Forecast_adjusted_scaled2) endif if (lobsdiag_forenkf) then deallocate(Observation_Operator_Jacobian_stind, Observation_Operator_Jacobian_endind, & Observation_Operator_Jacobian_val) endif enddo peloop ! ipe enddo ! satellite if (nanal == nanals) print *,'radiance ob profiles, total obs',nprof,nob if (nanal == nanals .and. lobsdiag_forenkf) print *,'time in calc_linhx for sat obs on proc',nproc,' = ',tsum if (nanal == nanals) print *,'time in read_raddiag_data for sat obs on proc',nproc,' = ',tsum2 if (nob /= nobs_max) then print *,'number of obs not what expected in get_satobs_data',nob,nobs_max call stop2(92) end if if (nobdiag /= nobs_maxdiag) then print *,'number of total diag obs not what expected in get_satobs_data',nobdiag,nobs_maxdiag call stop2(92) end if end subroutine get_satobs_data_nc ! write spread diagnostics subroutine write_satobs_data(obspath, datestring, nobs_max, nobs_maxdiag, x_fit, x_sprd, x_used, id, id2, gesid2) implicit none character*500, intent(in) :: obspath character(len=10), intent(in) :: datestring integer(i_kind), intent(in) :: nobs_max, nobs_maxdiag real(r_single), dimension(nobs_max), intent(in) :: x_fit, x_sprd integer(i_kind), dimension(nobs_maxdiag), intent(in) :: x_used character(len=10), intent(in) :: id, id2, gesid2 if (netcdf_diag) then call write_satobs_data_nc(obspath, datestring, nobs_max, nobs_maxdiag, x_fit, x_sprd, x_used, id, gesid2) else call write_satobs_data_bin(obspath, datestring, nobs_max, nobs_maxdiag, x_fit, x_sprd, x_used, id, id2, gesid2) endif end subroutine write_satobs_data ! write spread diagnostics to binary file subroutine write_satobs_data_bin(obspath, datestring, nobs_max, nobs_maxdiag, x_fit, x_sprd, x_used, id, id2, gesid2) use radinfo, only: iuse_rad,jpch_rad,nusis use read_diag, only: iversion_radiag_2, ireal_radiag, ireal_old_radiag implicit none character*500, intent(in) :: obspath character(len=10), intent(in) :: datestring integer(i_kind), intent(in) :: nobs_max, nobs_maxdiag real(r_single), dimension(nobs_max), intent(in) :: x_fit, x_sprd integer(i_kind), dimension(nobs_maxdiag), intent(in) :: x_used character(len=10), intent(in) :: id, id2, gesid2 character*500 obsfile,obsfile2 character(len=4) pe_name character(len=20) :: sat_type integer(i_kind) iunit,iunit2,iflag,nobs, nobsdiag,n,nsat,ipe,i,jpchstart logical fexist,init_pass character(len=10):: satid,sentype character(len=20):: sensat integer(i_kind):: jiter,nchanl,npred,ianldate,ireal,ipchan,iextra,jextra integer(i_kind):: idiag,angord,iversion,inewpc,isens,ijacob integer(i_kind):: iuse_tmp,nuchan_tmp,iochan_tmp real(r_single) :: freq_tmp,polar_tmp,wave_tmp,varch_tmp,tlapmean_tmp real(r_single),dimension(:,:),allocatable :: data_tmp real(r_single),dimension(:),allocatable :: fix_tmp real(r_single),dimension(:,:),allocatable :: extra_tmp iunit = 7 iunit2 = 17 nobs = 0 nobsdiag = 0 do nsat=1,nsats_rad jpchstart=0 do i=1,jpch_rad write(sat_type,'(a20)') adjustl(dsis(nsat)) ! The following is to sort out some historical naming conventions select case (sat_type(1:4)) case ('airs') sat_type='airs_aqua' case ('iasi') if (index(sat_type,'metop-a') /= 0) sat_type='iasi_metop-a' if (index(sat_type,'metop-b') /= 0) sat_type='iasi_metop-b' if (index(sat_type,'metop-c') /= 0) sat_type='iasi_metop-c' end select if(sat_type == trim(nusis(i)) .and. iuse_rad(i) > 0) then jpchstart = i exit end if end do if(jpchstart == 0) cycle init_pass = .true. if (datestring .eq. '0000000000') then obsfile2 = trim(adjustl(obspath))//"diag_"//trim(sattypes_rad(nsat))//"_"//trim(adjustl(gesid2))//"."//trim(adjustl(id2)) else obsfile2 = trim(adjustl(obspath))//"diag_"//trim(sattypes_rad(nsat))//"_"//trim(adjustl(gesid2))//"."//datestring//'_'//trim(adjustl(id2)) endif peloop: do ipe=0,npefiles write(pe_name,'(i4.4)') ipe if (npefiles .eq. 0) then ! diag file (concatenated pe* files) obsfile = trim(adjustl(obspath))//"diag_"//trim(sattypes_rad(nsat))//"_ges."//datestring//'_'//trim(adjustl(id)) inquire(file=obsfile,exist=fexist) if (.not. fexist .or. datestring .eq. '0000000000') then obsfile = trim(adjustl(obspath))//"diag_"//trim(sattypes_rad(nsat))//"_ges."//trim(adjustl(id)) endif else ! raw, unconcatenated pe* files. obsfile = trim(adjustl(obspath))//'gsitmp_'//trim(adjustl(id))//'/pe'//pe_name//'.'//trim(sattypes_rad(nsat))//'_01' endif inquire(file=obsfile,exist=fexist) if(.not.fexist) cycle peloop open(iunit,form="unformatted",file=obsfile) rewind(iunit) if (init_pass) then open(iunit2,form="unformatted",file=obsfile2) ! Read header (fixed_part). read(iunit,IOSTAT=iflag) sensat,satid,sentype,jiter,nchanl,npred,ianldate,& ireal,ipchan,iextra,jextra,idiag,angord,iversion,inewpc,isens,ijacob if (iflag == 0) then write(iunit2) sensat,satid,sentype,jiter,nchanl,npred,ianldate,& ireal,ipchan,iextra,jextra,idiag,angord,iversion,inewpc,isens,ijacob else rewind(iunit) read(iunit,IOSTAT=iflag) sensat,satid,sentype,jiter,nchanl,npred,ianldate,& ireal,ipchan,iextra,jextra,idiag,angord,iversion,inewpc,isens ijacob=0 if (iflag==0) then write(iunit2) sensat,satid,sentype,jiter,nchanl,npred,ianldate, & ireal,ipchan,iextra,jextra,idiag,angord,iversion,inewpc,isens else rewind(iunit) read(iunit,IOSTAT=iflag) sensat,satid,sentype,jiter,nchanl,npred,ianldate, & ireal,ipchan,iextra,jextra,idiag,angord,iversion,inewpc if (iflag==0) then write(iunit2) sensat,satid,sentype,jiter,nchanl,npred,ianldate, & ireal,ipchan,iextra,jextra,idiag,angord,iversion,inewpc else rewind(iunit) read(iunit,IOSTAT=iflag) sensat,satid,sentype,jiter,nchanl,npred,ianldate, & ireal,ipchan,iextra,jextra if (iflag==0) then write(iunit2) sensat,satid,sentype,jiter,nchanl,npred,ianldate, & ireal,ipchan,iextra,jextra else write(6,*)'READ_RADIAG_HEADER: ***ERROR*** Unknown file format.Cannot read' call stop2(5555) endif endif endif endif ! read header (channel part) do n=1, nchanl read(iunit,IOSTAT=iflag) freq_tmp,polar_tmp,wave_tmp,varch_tmp,tlapmean_tmp, & iuse_tmp,nuchan_tmp,iochan_tmp write(iunit2) freq_tmp,polar_tmp,wave_tmp,varch_tmp,tlapmean_tmp,iuse_tmp, & nuchan_tmp,iochan_tmp if (iflag/=0) return end do init_pass = .false. endif allocate(data_tmp(idiag,nchanl)) if (iversion < iversion_radiag_2) then allocate( fix_tmp( ireal_old_radiag ) ) else allocate( fix_tmp( ireal_radiag ) ) end if if (iextra > 0) then allocate(extra_tmp(iextra,jextra)) endif do if (iextra == 0) then read(iunit,IOSTAT=iflag) fix_tmp, data_tmp else read(iunit,IOSTAT=iflag) fix_tmp, data_tmp, extra_tmp endif if( iflag /= 0 ) then exit end if chan:do n=1,nchanl if (data_tmp(5,n) == 0) data_tmp(5,n) = 100 ! qcmark: not used in EnKF nobsdiag = nobsdiag + 1 if (x_used(nobsdiag) == 1) then nobs = nobs + 1 data_tmp(5,n) = 0 data_tmp(3,n) = x_fit(nobs) + data_tmp(3,n)-data_tmp(2,n) data_tmp(2,n) = x_fit(nobs) data_tmp(16+angord+3,n) = x_sprd(nobs) endif enddo chan if (iextra == 0) then write(iunit2) fix_tmp, data_tmp else write(iunit2) fix_tmp, data_tmp, extra_tmp endif enddo if (allocated(data_tmp)) deallocate(data_tmp) if (allocated(fix_tmp)) deallocate(fix_tmp) if (allocated(extra_tmp)) deallocate(extra_tmp) close(iunit) enddo peloop ! ipe close(iunit2) enddo ! satellite if (nobs /= nobs_max) then print *,'number of obs not what expected in get_satobs_data',nobs,nobs_max call stop2(92) end if if (nobsdiag /= nobs_maxdiag) then print *,'number of total diag obs not what expected in get_satobs_data',nobsdiag,nobs_maxdiag call stop2(92) end if end subroutine write_satobs_data_bin ! writing spread diagnostics to netcdf file subroutine write_satobs_data_nc(obspath, datestring, nobs_max, nobs_maxdiag, & x_fit, x_sprd, x_used, id, gesid) use netcdf, only: nf90_inq_dimid, nf90_open, nf90_close, NF90_NETCDF4, & nf90_inquire_dimension, NF90_WRITE, NF90_NOWRITE, nf90_create, nf90_def_dim use ncdw_climsg, only: nclayer_check use radinfo, only: iuse_rad,nusis,jpch_rad use constants, only: r_missing implicit none character*500, intent(in) :: obspath character(len=10), intent(in) :: datestring integer(i_kind), intent(in) :: nobs_max, nobs_maxdiag real(r_single), dimension(nobs_max), intent(in) :: x_fit, x_sprd integer(i_kind), dimension(nobs_maxdiag), intent(in) :: x_used character(len=10), intent(in) :: id, gesid character*500 obsfile, obsfile2 character(len=4) pe_name character(len=20) :: sat_type integer(i_kind) :: iunit, nobsid integer(i_kind) :: nob, nobdiag, nobs, ipe, i, nsat, jpchstart integer(i_kind), dimension(:), allocatable :: enkf_use_flag real(r_single), dimension(:), allocatable :: enkf_fit, enkf_sprd logical :: fexist nob = 0 nobdiag = 0 do nsat=1,nsats_rad jpchstart=0 do i=1,jpch_rad write(sat_type,'(a20)') adjustl(dsis(nsat)) ! The following is to sort out some historical naming conventions select case (sat_type(1:4)) case ('airs') sat_type='airs_aqua' case ('iasi') if (index(sat_type,'metop-a') /= 0) sat_type='iasi_metop-a' if (index(sat_type,'metop-b') /= 0) sat_type='iasi_metop-b' if (index(sat_type,'metop-c') /= 0) sat_type='iasi_metop-c' end select if(sat_type == trim(nusis(i)) .and. iuse_rad(i) > 0) then jpchstart = i exit end if end do if(jpchstart == 0) cycle peloop: do ipe=0,npefiles write(pe_name,'(i4.4)') ipe if (npefiles .eq. 0) then ! diag file (concatenated pe* files) obsfile = trim(adjustl(obspath))//"diag_"//trim(sattypes_rad(nsat))//"_ges."//datestring//'_'//trim(adjustl(id))//".nc4" obsfile2 = trim(adjustl(obspath))//"diag_"//trim(sattypes_rad(nsat))//"_ges."//datestring//'_'//trim(adjustl(id))//"_spread.nc4" inquire(file=obsfile,exist=fexist) if (.not. fexist .or. datestring .eq. '0000000000') then obsfile = trim(adjustl(obspath))//"diag_"//trim(sattypes_rad(nsat))//"_ges."//trim(adjustl(id))//".nc4" obsfile2 = trim(adjustl(obspath))//"diag_"//trim(sattypes_rad(nsat))//"_ges."//trim(adjustl(id))//"_spread.nc4" endif else ! raw, unconcatenated pe* files. obsfile = trim(adjustl(obspath))//'gsitmp_'//trim(adjustl(id))//'/pe'//pe_name//'.'//trim(sattypes_rad(nsat))//'_01'//".nc4" obsfile2 = trim(adjustl(obspath))//'gsitmp_'//trim(adjustl(id))//'/pe'//pe_name//'.'//trim(sattypes_rad(nsat))//'_01'//"_spread.nc4" endif inquire(file=obsfile,exist=fexist) if (.not. fexist) cycle peloop call nclayer_check(nf90_open(obsfile, NF90_NOWRITE, iunit)) call nclayer_check(nf90_inq_dimid(iunit, "nobs", nobsid)) call nclayer_check(nf90_inquire_dimension(iunit, nobsid, len = nobs)) call nclayer_check(nf90_close(iunit)) if (nobs <= 0) cycle peloop allocate(enkf_use_flag(nobs), enkf_fit(nobs), enkf_sprd(nobs)) do i = 1, nobs nobdiag = nobdiag + 1 ! skip if not used in EnKF if (x_used(nobdiag) == 1) then ! update if it is used in EnKF nob = nob + 1 enkf_use_flag(i) = 1 enkf_fit(i) = x_fit(nob) enkf_sprd(i) = x_sprd(nob) else enkf_use_flag(i) = -1 enkf_fit(i) = r_missing enkf_sprd(i) = r_missing endif enddo inquire(file=obsfile2,exist=fexist) if (.not. fexist) then call nclayer_check(nf90_create(trim(obsfile2), NF90_NETCDF4, & iunit)) call nclayer_check(nf90_def_dim(iunit, "nobs", nobs, nobsid)) else call nclayer_check(nf90_open(obsfile2, NF90_WRITE, iunit)) call nclayer_check(nf90_inq_dimid(iunit, "nobs", nobsid)) endif call write_ncvar_int(iunit, nobsid, "EnKF_use_flag", enkf_use_flag) call write_ncvar_single(iunit, nobsid, "EnKF_fit_"//trim(gesid), enkf_fit) call write_ncvar_single(iunit, nobsid, "EnKF_spread_"//trim(gesid), enkf_sprd) call nclayer_check(nf90_close(iunit)) deallocate(enkf_use_flag, enkf_fit, enkf_sprd) enddo peloop ! ipe loop enddo if (nob .ne. nobs_max) then print *,'number of obs not what expected in write_satobs_data',nob,nobs_max call stop2(94) end if if (nobdiag /= nobs_maxdiag) then print *,'number of total obs in diag not what expected in write_satobs_data',nobdiag, nobs_maxdiag call stop2(94) endif contains subroutine write_ncvar_single(iunit, dimid, varname, field) use netcdf, only: nf90_def_var, nf90_put_var, nf90_inq_varid, & nf90_def_var_deflate,NF90_FLOAT, NF90_ENOTVAR use ncdw_climsg, only: nclayer_check use ncdw_types, only: NLAYER_COMPRESSION implicit none integer(i_kind), intent(in) :: iunit, dimid character(*), intent(in) :: varname real(r_single), dimension(:), allocatable :: field integer :: ierr, varid ierr = nf90_inq_varid(iunit, varname, varid) if (ierr == NF90_ENOTVAR) then call nclayer_check(nf90_def_var(iunit, varname, NF90_FLOAT, dimid, varid)) call nclayer_check(nf90_def_var_deflate(iunit, varid, 1, 1, int(NLAYER_COMPRESSION))) endif call nclayer_check(nf90_put_var(iunit, varid, field)) end subroutine write_ncvar_single subroutine write_ncvar_int(iunit, dimid, varname, field) use netcdf, only: nf90_def_var, nf90_put_var, nf90_inq_varid, & nf90_def_var_deflate,NF90_INT, NF90_ENOTVAR use ncdw_climsg, only: nclayer_check use ncdw_types, only: NLAYER_COMPRESSION implicit none integer(i_kind), intent(in) :: iunit, dimid character(*), intent(in) :: varname integer(i_kind), dimension(:), allocatable :: field integer :: ierr, varid ierr = nf90_inq_varid(iunit, varname, varid) if (ierr == NF90_ENOTVAR) then call nclayer_check(nf90_def_var(iunit, varname, NF90_INT, dimid, varid)) call nclayer_check(nf90_def_var_deflate(iunit, varid, 1, 1, int(NLAYER_COMPRESSION))) endif call nclayer_check(nf90_put_var(iunit, varid, field)) end subroutine write_ncvar_int end subroutine write_satobs_data_nc end module readsatobs