From de74bb2a5a14b61b6131a242c8bc600e1936aeb9 Mon Sep 17 00:00:00 2001 From: James Jung Date: Tue, 30 Jul 2024 20:16:10 -0400 Subject: [PATCH] Format changes for EUMETSAT metop-sg and CADS debug fix (#773) --- src/gsi/gsimod.F90 | 2 +- src/gsi/mrmsmod.f90 | 7 +++- src/gsi/obsmod.F90 | 3 +- src/gsi/qcmod.f90 | 2 +- src/gsi/radinfo.f90 | 4 +- src/gsi/read_cris.f90 | 89 +++++++++++++++++++++++------------------- src/gsi/read_diag.f90 | 8 ++-- src/gsi/read_iasi.f90 | 81 +++++++++++++++++++++----------------- src/gsi/read_obs.F90 | 59 ++++++++++++++-------------- src/gsi/read_viirs.f90 | 2 +- src/gsi/statsrad.f90 | 6 +-- 11 files changed, 144 insertions(+), 119 deletions(-) diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index d7f5667252..8c6977d78c 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -2320,7 +2320,7 @@ subroutine gsimain_initialize endif do i=1,ndat write(6,401)dfile(i),dtype(i),dplat(i),dsis(i),dval(i),dthin(i),dsfcalc(i),time_window(i) - 401 format(1x,a20,1x,a10,1x,a10,1x,a20,1x,f10.2,1x,I3,1x,I3,1x,f10.2) + 401 format(1x,a20,1x,a10,1x,a11,1x,a20,1x,f10.2,1x,I3,1x,I3,1x,f10.2) end do write(6,superob_radar) write(6,lag_data) diff --git a/src/gsi/mrmsmod.f90 b/src/gsi/mrmsmod.f90 index 713b974d7a..1cbad54904 100644 --- a/src/gsi/mrmsmod.f90 +++ b/src/gsi/mrmsmod.f90 @@ -27,7 +27,9 @@ subroutine load_mrms_data_info (mrms_listfile,nrows0,ntot_mrms,nrows_mrms,nrows, integer(i_kind),parameter::nobstype_mrms=2 ! first for vr, and the second for ref integer(i_kind),intent(in):: nrows0,ntot_mrms,nrows_mrms,nrows character(len=*),intent(in),optional ::mrms_listfile, rcname ! input filename - character(10),intent(inout),dimension(nrows):: dtype,ditype,dplat + character(len=*),intent(inout),dimension(nrows):: dtype + character(len=*),intent(inout),dimension(nrows):: ditype + character(len=*),intent(inout),dimension(nrows):: dplat character(20),intent(inout),dimension(nrows):: obsfile_all character(*),intent(inout),dimension(nrows):: dfile character(20),intent(inout),dimension(nrows):: dsis @@ -42,7 +44,8 @@ subroutine load_mrms_data_info (mrms_listfile,nrows0,ntot_mrms,nrows_mrms,nrows, real(r_kind),allocatable,dimension(:):: dmesh_mrms - character(10),allocatable,dimension(:):: dtype_mrms,ditype_mrms,dplat_mrms + character(10),allocatable,dimension(:):: dtype_mrms,ditype_mrms + character(11),allocatable,dimension(:):: dplat_mrms character(120),allocatable,dimension(:):: dfile_mrms character(20),allocatable,dimension(:):: dsis_mrms real(r_kind) ,allocatable,dimension(:):: dval_mrms diff --git a/src/gsi/obsmod.F90 b/src/gsi/obsmod.F90 index 19f6a072fb..eab9c828f7 100644 --- a/src/gsi/obsmod.F90 +++ b/src/gsi/obsmod.F90 @@ -626,7 +626,8 @@ module obsmod character(128) dirname character(128) obs_input_common character(20),allocatable,dimension(:):: obsfile_all - character(10),allocatable,dimension(:):: dtype,ditype,dplat + character(10),allocatable,dimension(:):: dtype,ditype + character(11),allocatable,dimension(:):: dplat character(120),allocatable,dimension(:):: dfile character(20),allocatable,dimension(:):: dsis real(r_kind) ,allocatable,dimension(:):: dval diff --git a/src/gsi/qcmod.f90 b/src/gsi/qcmod.f90 index 9804965573..2c4cbf6317 100644 --- a/src/gsi/qcmod.f90 +++ b/src/gsi/qcmod.f90 @@ -4240,7 +4240,7 @@ subroutine qc_goesimg(nchanl,is,ndat,nsig,ich,dplat,sea,land,ice,snow,luse, & real(r_kind),dimension(nsig,nchanl),intent(in ) :: temp,wmix real(r_kind),dimension(nchanl), intent(in ) :: tb_obs,tb_obs_sdv,tbc,tnoise,emissivity_k,ts real(r_kind),dimension(nchanl), intent(inout) :: errf,varinv - character(10), intent(in ) :: dplat + character(len=*), intent(in ) :: dplat ! Declare local parameters diff --git a/src/gsi/radinfo.f90 b/src/gsi/radinfo.f90 index 4ad17626e6..8bb496c645 100644 --- a/src/gsi/radinfo.f90 +++ b/src/gsi/radinfo.f90 @@ -805,10 +805,10 @@ subroutine radinfo_read end do close(lunin) 100 format(a1,a120) -110 format(i4,1x,a20,' chan= ',i4, & +110 format(i4,1x,a20,' chan= ',i5, & ' var= ',f7.3,' varch_cld=',f7.3,' use= ',i2,' ermax= ',F7.3, & ' b_rad= ',F7.2,' pg_rad=',F7.2,' icld_det=',I2,' icloud=',I2,' iaeros=',I2) -111 format(i4,1x,a20,' chan= ',i4, & +111 format(i4,1x,a20,' chan= ',i5, & ' var= ',f7.3,' varch_cld=',f7.3,' use= ',i2,' ermax= ',F7.3, & ' b_rad= ',F7.2,' pg_rad=',F7.2,' iextra_det=',I2, 'icloud=',I2,'iaeros=', I2) diff --git a/src/gsi/read_cris.f90 b/src/gsi/read_cris.f90 index d5668f8864..84288a7f04 100644 --- a/src/gsi/read_cris.f90 +++ b/src/gsi/read_cris.f90 @@ -112,9 +112,9 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& integer(i_kind) ,intent(in ) :: mype_sub integer(i_kind) ,intent(in ) :: npe_sub integer(i_kind) ,intent(in ) :: mpi_comm_sub - character(len=*), intent(in ) :: infile - character(len=10),intent(in ) :: jsatid - character(len=*), intent(in ) :: obstype + character(len=*) ,intent(in ) :: infile + character(len=*) ,intent(in ) :: jsatid + character(len=*) ,intent(in ) :: obstype character(len=20),intent(in ) :: sis real(r_kind) ,intent(in ) :: twind real(r_kind) ,intent(inout) :: val_cris @@ -200,7 +200,6 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& real(r_kind),dimension(83,7) :: imager_info real(r_kind),dimension(7) :: imager_cluster_size real(r_kind),dimension(2) :: imager_mean, imager_std_dev, imager_conversion - real(r_kind) :: imager_cluster_tot ! bufr error codes ! real(r_kind),dimension(7,3) :: error_codes @@ -848,10 +847,9 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& if ( cris_cads ) then call ufbseq(lnbufr,imager_info,83,7,iret,'CRISCS') if ( iret == 7 .and. imager_info(3,1) <= 100.0_r_kind .and. & - imager_info(3,1) >= zero .and. imager_coeff ) then ! if imager cluster info exists + sum(imager_info(3,:)) > zero .and. imager_coeff ) then ! if imager cluster info exists imager_mean = zero imager_std_dev = zero - imager_cluster_tot = zero imager_cluster_flag = .TRUE. imager_cluster_size = imager_info(3,1:7) imager_cluster_size(:) = imager_cluster_size(:) / sum(imager_cluster_size(:)) @@ -871,51 +869,62 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& imager_cluster_info: do j=1,7 i = imager_cluster_index(j) - data_all(maxinfo+j,itx) = imager_cluster_size(i) ! Imager cluster fraction - imager_cluster_tot = imager_cluster_tot + imager_info(3,i) +! If the cluster size, or radiance values of channel 4 or 5 are zero, do not compute statistics for that cluster - iexponent = -(nint(imager_info(75,i)) -11) ! channel 15 radiance for each cluster - imager_info(76,i) = imager_info(76,i) * imager_conversion(1) * (ten ** iexponent) + if ( imager_cluster_size(i) > zero .and. imager_info(76,i) > zero .and. imager_info(81,i) > zero ) then + data_all(maxinfo+j,itx) = imager_cluster_size(i) ! Imager cluster fraction - iexponent = -(nint(imager_info(77,i)) -11) ! channel 15 radiance std dev for each cluster. - imager_info(78,i) = imager_info(78,i) * imager_conversion(1) * (ten ** iexponent) + iexponent = -(nint(imager_info(75,i)) -11) ! channel 15 radiance for each cluster + imager_info(76,i) = imager_info(76,i) * imager_conversion(1) * (ten ** iexponent) - call crtm_planck_temperature(sensorindex_imager,4,imager_info(76,i),data_all(maxinfo+7+j,itx)) - data_all(maxinfo+7+j,itx) = max(data_all(maxinfo+7+j,itx),zero) + iexponent = -(nint(imager_info(77,i)) -11) ! channel 15 radiance std dev for each cluster. + imager_info(78,i) = imager_info(78,i) * imager_conversion(1) * (ten ** iexponent) - iexponent = -(nint(imager_info(80,i)) -11) ! channel 16 radiance for each cluster - imager_info(81,i) = imager_info(81,i) * imager_conversion(2) * (ten ** iexponent) + iexponent = -(nint(imager_info(80,i)) -11) ! channel 16 radiance for each cluster + imager_info(81,i) = imager_info(81,i) * imager_conversion(2) * (ten ** iexponent) - iexponent = -(nint(imager_info(82,i))-5 ) ! channel 16 radiance std dev for each cluster. - iexponent = -(nint(imager_info(82,i)) -11) ! channel 16 radiance std dev for each cluster. - imager_info(83,i) = imager_info(83,i) * imager_conversion(2) * (ten ** iexponent) + iexponent = -(nint(imager_info(82,i)) -11) ! channel 16 radiance std dev for each cluster. + imager_info(83,i) = imager_info(83,i) * imager_conversion(2) * (ten ** iexponent) - call crtm_planck_temperature(sensorindex_imager,5,imager_info(81,i),data_all(maxinfo+14+j,itx)) - data_all(maxinfo+14+j,itx) = max(data_all(maxinfo+14+j,itx),zero) - - - end do imager_cluster_info + call crtm_planck_temperature(sensorindex_imager,4,imager_info(76,i),data_all(maxinfo+7+j,itx)) + data_all(maxinfo+7+j,itx) = max(data_all(maxinfo+7+j,itx),zero) + call crtm_planck_temperature(sensorindex_imager,5,imager_info(81,i),data_all(maxinfo+14+j,itx)) + data_all(maxinfo+14+j,itx) = max(data_all(maxinfo+14+j,itx),zero) + else + data_all(maxinfo+j,itx) = zero ! something is wrong + data_all(maxinfo+7+j,itx) = zero ! set everything to zero + data_all(maxinfo+14+j,itx) = zero + endif + end do imager_cluster_info ! Compute cluster averages for each channel - imager_mean(1) = sum(imager_cluster_size(:) * imager_info(76,:)) ! Channel 15 radiance cluster average - imager_std_dev(1) = sum(imager_cluster_size(:) * (imager_info(76,:)**2 + imager_info(78,:)**2)) - imager_mean(1)**2 - imager_std_dev(1) = sqrt(max(imager_std_dev(1),zero)) ! Channel 15 radiance RMSE - call crtm_planck_temperature(sensorindex_imager,4,(imager_std_dev(1) + imager_mean(1)),imager_std_dev(1)) - call crtm_planck_temperature(sensorindex_imager,4,imager_mean(1),imager_mean(1)) ! Channel 15 average BT - imager_std_dev(1) = imager_std_dev(1) - imager_mean(1) ! Channel 15 BT std dev - data_all(maxinfo+22,itx) = imager_std_dev(1) - - imager_mean(2) = sum(imager_cluster_size(:) * imager_info(81,:)) ! Channel 16 radiance cluster average - imager_std_dev(2) = sum(imager_cluster_size(:) * (imager_info(81,:)**2 + imager_info(83,:)**2)) - imager_mean(1)**2 - imager_std_dev(2) = sqrt(max(imager_std_dev(1),zero)) ! Channel 16 radiance RMSE - call crtm_planck_temperature(sensorindex_imager,5,(imager_std_dev(2) + imager_mean(2)),imager_std_dev(2)) - call crtm_planck_temperature(sensorindex_imager,5,imager_mean(2),imager_mean(2)) ! Channel 16 average BT - imager_std_dev(2) = imager_std_dev(2) - imager_mean(2) ! Channel 16 BT std dev - data_all(maxinfo+23,itx) = imager_std_dev(2) + imager_mean(1) = sum(imager_cluster_size(:) * imager_info(76,:)) ! Channel 15 radiance cluster average + imager_std_dev(1) = sum(imager_cluster_size(:) * (imager_info(76,:)**2 + imager_info(78,:)**2)) - imager_mean(1)**2 + imager_std_dev(1) = sqrt(max(imager_std_dev(1),zero)) ! Channel 15 radiance RMSE + if ( imager_mean(1) > zero .and. imager_std_dev(1) > zero ) then + call crtm_planck_temperature(sensorindex_imager,4,(imager_std_dev(1) + imager_mean(1)),imager_std_dev(1)) + call crtm_planck_temperature(sensorindex_imager,4,imager_mean(1),imager_mean(1)) ! Channel 15 average BT + imager_std_dev(1) = imager_std_dev(1) - imager_mean(1) ! Channel 15 BT std dev + data_all(maxinfo+22,itx) = imager_std_dev(1) + else + data_all(maxinfo+22,itx) = zero + endif + + imager_mean(2) = sum(imager_cluster_size(:) * imager_info(81,:)) ! Channel 16 radiance cluster average + imager_std_dev(2) = sum(imager_cluster_size(:) * (imager_info(81,:)**2 + imager_info(83,:)**2)) - imager_mean(1)**2 + imager_std_dev(2) = sqrt(max(imager_std_dev(1),zero)) ! Channel 16 radiance RMSE + if ( imager_mean(2) > zero .and. imager_std_dev(2) > zero ) then + call crtm_planck_temperature(sensorindex_imager,5,(imager_std_dev(2) + imager_mean(2)),imager_std_dev(2)) + call crtm_planck_temperature(sensorindex_imager,5,imager_mean(2),imager_mean(2)) ! Channel 16 average BT + imager_std_dev(2) = imager_std_dev(2) - imager_mean(2) ! Channel 16 BT std dev + data_all(maxinfo+23,itx) = imager_std_dev(2) + else + data_all(maxinfo+23,itx) = zero + endif else ! Imager cluster information is missing. Set everything to zero - data_all(maxinfo+1 : maxinfo+25,itx) = zero + data_all(maxinfo+1 : maxinfo+cads_info,itx) = zero endif endif ! cris_cads diff --git a/src/gsi/read_diag.f90 b/src/gsi/read_diag.f90 index e389a708d5..f69f2f5133 100644 --- a/src/gsi/read_diag.f90 +++ b/src/gsi/read_diag.f90 @@ -88,7 +88,7 @@ module read_diag ! Declare structures for radiance diagnostic file information type diag_header_fix_list character(len=20) :: isis ! sat and sensor type - character(len=10) :: id ! sat type + character(len=11) :: id ! sat type character(len=10) :: obstype ! observation type integer(i_kind) :: jiter ! outer loop counter integer(i_kind) :: nchan ! number of channels in the sensor @@ -414,7 +414,8 @@ subroutine read_radiag_header_nc(ftin,header_fix,header_chan,iflag) real(r_kind),allocatable,dimension(:) :: r_var_stor integer(i_kind),allocatable,dimension(:) :: i_var_stor character(20) :: isis - character(10) :: id, obstype + character(11) :: id + character(10) :: obstype ! integer(i_kind),dimension(:),allocatable :: jiter, nchan_diag, npred, idate, & integer(i_kind) :: jiter, nchan_diag, npred, idate, & ireal, ipchan, iextra, jextra, & @@ -515,7 +516,8 @@ subroutine read_radiag_header_bin(ftin,npred_radiag,retrieval,header_fix,header_ ! Declare local variables character(len=2):: string - character(len=10):: satid,sentype + character(len=11):: satid + character(len=10):: sentype character(len=20):: sensat integer(i_kind) :: i,ich integer(i_kind):: jiter,nchanl,npred,ianldate,ireal,ipchan,iextra,jextra diff --git a/src/gsi/read_iasi.f90 b/src/gsi/read_iasi.f90 index edd7a9b50e..75dc50bb76 100644 --- a/src/gsi/read_iasi.f90 +++ b/src/gsi/read_iasi.f90 @@ -142,9 +142,9 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& integer(i_kind) ,intent(in ) :: mype_sub integer(i_kind) ,intent(in ) :: npe_sub integer(i_kind) ,intent(in ) :: mpi_comm_sub - character(len=*), intent(in ) :: infile - character(len=10),intent(in ) :: jsatid - character(len=*), intent(in ) :: obstype + character(len=*) ,intent(in ) :: infile + character(len=*) ,intent(in ) :: jsatid + character(len=*) ,intent(in ) :: obstype character(len=20),intent(in ) :: sis real(r_kind) ,intent(in ) :: twind real(r_kind) ,intent(inout) :: val_iasi @@ -233,7 +233,6 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& real(r_kind),dimension(33,7) :: imager_info real(r_kind),dimension(7) :: imager_cluster_size real(r_kind),dimension(2) :: imager_mean, imager_std_dev - real(r_kind) :: imager_cluster_tot ! Set standard parameters character(8),parameter:: fov_flag="crosstrk" @@ -813,10 +812,9 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& if ( iasi_cads ) then call ufbseq(lnbufr,imager_info,33,7,iret,'IASIL1CS') if (iret == 7 .and. imager_info(3,1) <= 100.0_r_kind .and. & - imager_info(3,1) >= zero .and. imager_coeff ) then ! if imager cluster info exists + sum(imager_info(3,:)) > zero .and. imager_coeff ) then ! if imager cluster info exists imager_mean = zero imager_std_dev = zero - imager_cluster_tot = zero imager_cluster_flag = .TRUE. imager_cluster_size = imager_info(3,1:7) imager_cluster_size(:) = imager_cluster_size(:) / sum(imager_cluster_size(:)) @@ -834,49 +832,62 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& imager_cluster_info: do j=1,7 i = imager_cluster_index(j) - data_all(maxinfo+j,itx) = imager_cluster_size(i) ! Imager cluster fraction - imager_cluster_tot = imager_cluster_tot + imager_info(3,i) +! If the cluster size, or radiance values of channel 4 or 5 are zero, do not compute statistics for that cluster + if ( imager_cluster_size(i) > zero .and. imager_info(26,i) > zero .and. imager_info(31,i) > zero ) then + data_all(maxinfo+j,itx) = imager_cluster_size(i) ! Imager cluster fraction - iexponent = -(nint(imager_info(25,i))-5 ) ! channel 4 radiance for each cluster. - imager_info(26,i) = imager_info(26,i) * (ten ** iexponent) + iexponent = -(nint(imager_info(25,i))-5 ) ! channel 4 radiance for each cluster. + imager_info(26,i) = imager_info(26,i) * (ten ** iexponent) - iexponent = -(nint(imager_info(27,i))-5 ) ! channel 4 radiance std dev for each cluster. - imager_info(28,i) = imager_info(28,i) * (ten ** iexponent) + iexponent = -(nint(imager_info(27,i))-5 ) ! channel 4 radiance std dev for each cluster. + imager_info(28,i) = imager_info(28,i) * (ten ** iexponent) - call crtm_planck_temperature(sensorindex_imager,2,imager_info(26,i),data_all(maxinfo+7+j,itx)) - data_all(maxinfo+7+j,itx) = max(data_all(maxinfo+7+j,itx),zero) + iexponent = -(nint(imager_info(30,i))-5 ) ! channel 5 radiance for each cluster + imager_info(31,i) = imager_info(31,i) * (ten ** iexponent) - iexponent = -(nint(imager_info(30,i))-5 ) ! channel 5 radiance for each cluster - imager_info(31,i) = imager_info(31,i) * (ten ** iexponent) + iexponent = -(nint(imager_info(32,i))-5 ) ! channel 5 radiance std dev for each cluser. + imager_info(33,i) = imager_info(33,i) * (ten ** iexponent) - iexponent = -(nint(imager_info(32,i))-5 ) ! channel 5 radiance std dev for each cluser. - imager_info(33,i) = imager_info(33,i) * (ten ** iexponent) - - call crtm_planck_temperature(sensorindex_imager,3,imager_info(31,i),data_all(maxinfo+14+j,itx)) - data_all(maxinfo+14+j,itx) = max(data_all(maxinfo+14+j,itx),zero) + call crtm_planck_temperature(sensorindex_imager,2,imager_info(26,i),data_all(maxinfo+7+j,itx)) + data_all(maxinfo+7+j,itx) = max(data_all(maxinfo+7+j,itx),zero) + call crtm_planck_temperature(sensorindex_imager,3,imager_info(31,i),data_all(maxinfo+14+j,itx)) + data_all(maxinfo+14+j,itx) = max(data_all(maxinfo+14+j,itx),zero) + else ! something is wrong + data_all(maxinfo+j,itx) = zero ! set everything to zero + data_all(maxinfo+7+j,itx) = zero + data_all(maxinfo+14+j,itx) = zero + endif end do imager_cluster_info ! Compute cluster averages for each channel - imager_mean(1) = sum(imager_cluster_size(:) * imager_info(26,:)) ! Channel 4 radiance cluster average + imager_mean(1) = sum(imager_cluster_size(:) * imager_info(26,:)) ! Channel 4 radiance cluster average imager_std_dev(1) = sum(imager_cluster_size(:) * (imager_info(26,:)**2 + imager_info(28,:)**2)) - imager_mean(1)**2 - imager_std_dev(1) = sqrt(max(imager_std_dev(1),zero)) ! Channel 4 radiance RMSE - call crtm_planck_temperature(sensorindex_imager,2,(imager_std_dev(1) + imager_mean(1)),imager_std_dev(1)) - call crtm_planck_temperature(sensorindex_imager,2,imager_mean(1),imager_mean(1)) ! Channel 4 average BT - imager_std_dev(1) = imager_std_dev(1) - imager_mean(1) ! Channel 4 BT std dev - data_all(maxinfo+22,itx) = imager_std_dev(1) - - imager_mean(2) = sum(imager_cluster_size(:) * imager_info(31,:)) ! Channel 5 radiance cluster average + imager_std_dev(1) = sqrt(max(imager_std_dev(1),zero)) ! Channel 4 radiance RMSE + if ( imager_mean(1) > zero .and. imager_std_dev(1) > zero ) then + call crtm_planck_temperature(sensorindex_imager,2,(imager_std_dev(1) + imager_mean(1)),imager_std_dev(1)) + call crtm_planck_temperature(sensorindex_imager,2,imager_mean(1),imager_mean(1)) ! Channel 4 average BT + imager_std_dev(1) = imager_std_dev(1) - imager_mean(1) ! Channel 4 BT std dev + data_all(maxinfo+22,itx) = imager_std_dev(1) + else + data_all(maxinfo+22,itx) = zero + endif + + imager_mean(2) = sum(imager_cluster_size(:) * imager_info(31,:)) ! Channel 5 radiance cluster average imager_std_dev(2) = sum(imager_cluster_size(:) * (imager_info(31,:)**2 + imager_info(33,:)**2)) - imager_mean(1)**2 - imager_std_dev(2) = sqrt(max(imager_std_dev(1),zero)) ! Channel 5 radiance RMSE - call crtm_planck_temperature(sensorindex_imager,3,(imager_std_dev(2) + imager_mean(2)),imager_std_dev(2)) - call crtm_planck_temperature(sensorindex_imager,3,imager_mean(2),imager_mean(2)) ! Channel 5 average BT - imager_std_dev(2) = imager_std_dev(2) - imager_mean(2) ! Channel 5 BT std dev - data_all(maxinfo+23,itx) = imager_std_dev(2) + imager_std_dev(2) = sqrt(max(imager_std_dev(1),zero)) ! Channel 5 radiance RMSE + if ( imager_mean(2) > zero .and. imager_std_dev(2) > zero ) then + call crtm_planck_temperature(sensorindex_imager,3,(imager_std_dev(2) + imager_mean(2)),imager_std_dev(2)) + call crtm_planck_temperature(sensorindex_imager,3,imager_mean(2),imager_mean(2)) ! Channel 5 average BT + imager_std_dev(2) = imager_std_dev(2) - imager_mean(2) ! Channel 5 BT std dev + data_all(maxinfo+23,itx) = imager_std_dev(2) + else + data_all(maxinfo+23,itx) = zero + endif else ! Imager cluster information is missing. Set everything to zero - data_all(maxinfo+1 : maxinfo+25,itx) = zero + data_all(maxinfo+1 : maxinfo+cads_info,itx) = zero endif endif ! iasi_cads = .true. ! diff --git a/src/gsi/read_obs.F90 b/src/gsi/read_obs.F90 index 6f0e32a4b6..b81bdb5752 100644 --- a/src/gsi/read_obs.F90 +++ b/src/gsi/read_obs.F90 @@ -791,7 +791,7 @@ subroutine read_obs(ndata,mype) logical,dimension(ndat):: belong,parallel_read,ears_possible,db_possible logical :: modis,use_sfc_any logical :: acft_profl_file - character(10):: obstype,platid + character(10):: obstype character(22):: string character(120):: infile character(20):: sis @@ -1415,7 +1415,6 @@ subroutine read_obs(ndata,mype) task_belongs: & if (i > 0 .and. belong(i)) then - platid=dplat(i) ! platid - satellites to read obstype=dtype(i) ! obstype - observation types to process infile=trim(dfile(i)) ! infile - units from which to read data sis=dsis(i) ! sensor/instrument/satellite indicator @@ -1553,12 +1552,12 @@ subroutine read_obs(ndata,mype) ! Process conventional SST (nsstbufr, at this moment) data elseif ( obstype == 'sst' ) then - string="--"//trim(ditype(i))//":sst:"//trim(platid) - if ( platid == 'nsst') then + string="--"//trim(ditype(i))//":sst:"//trim(dplat(i)) + if ( dplat(i) == 'nsst') then call read_nsstbufr(nread,npuse,nouse,gstime,infile,obstype, & lunout,twind,sis,nobs_sub1(1,i)) string='READ_NSSTBUFR' - elseif ( platid == 'mods') then + elseif ( dplat(i) == 'mods') then call read_modsbufr(nread,npuse,nouse,gstime,infile,obstype, & lunout,twind,sis,nobs_sub1(1,i)) string='READ_MODSBUFR' @@ -1716,12 +1715,12 @@ subroutine read_obs(ndata,mype) ! Process TOVS 1b data rad_obstype_select: & - if (platid /= 'aqua' .and. (obstype == 'amsua' .or. & + if (dplat(i) /= 'aqua' .and. (obstype == 'amsua' .or. & obstype == 'amsub' .or. obstype == 'msu' .or. & obstype == 'mhs' .or. obstype == 'hirs4' .or. & obstype == 'hirs3' .or. obstype == 'hirs2' .or. & obstype == 'ssu' )) then - call read_bufrtovs(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,& + call read_bufrtovs(mype,val_dat,ithin,isfcalc,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis, & mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), nobs_sub1(1,i), & read_rec(i),read_ears_rec(i),read_db_rec(i),dval_use,radmod) @@ -1729,7 +1728,7 @@ subroutine read_obs(ndata,mype) ! Process atms data else if (obstype == 'atms') then - call read_atms(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,& + call read_atms(mype,val_dat,ithin,isfcalc,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis, & mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i),nobs_sub1(1,i),& read_rec(i),read_ears_rec(i),read_db_rec(i),dval_use,radmod) @@ -1737,7 +1736,7 @@ subroutine read_obs(ndata,mype) ! Process saphir data else if (obstype == 'saphir') then - call read_saphir(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,& + call read_saphir(mype,val_dat,ithin,isfcalc,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis, & mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i),dval_use) @@ -1745,9 +1744,9 @@ subroutine read_obs(ndata,mype) ! Process airs data - else if(platid == 'aqua' .and. (obstype == 'airs' .or. & + else if(dplat(i) == 'aqua' .and. (obstype == 'airs' .or. & obstype == 'amsua' .or. obstype == 'hsb' ))then - call read_airs(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,& + call read_airs(mype,val_dat,ithin,isfcalc,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis,& mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i),read_rec(i),dval_use) @@ -1755,7 +1754,7 @@ subroutine read_obs(ndata,mype) ! Process iasi data else if(obstype == 'iasi')then - call read_iasi(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,& + call read_iasi(mype,val_dat,ithin,isfcalc,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis,& mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i),nobs_sub1(1,i), & read_rec(i),read_ears_rec(i),read_db_rec(i),dval_use) @@ -1763,7 +1762,7 @@ subroutine read_obs(ndata,mype) ! Process cris data else if(obstype == 'cris' .or. obstype =='cris-fsr' )then - call read_cris(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,& + call read_cris(mype,val_dat,ithin,isfcalc,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis,& mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i),nobs_sub1(1,i), & read_rec(i),read_ears_rec(i),read_db_rec(i),dval_use) @@ -1774,7 +1773,7 @@ subroutine read_obs(ndata,mype) else if (obstype == 'sndr' .or. & obstype == 'sndrd1' .or. obstype == 'sndrd2' .or. & obstype == 'sndrd3' .or. obstype == 'sndrd4') then - call read_goesndr(mype,val_dat,ithin,rmesh,platid,& + call read_goesndr(mype,val_dat,ithin,rmesh,dplat(i),& infile,lunout,obstype,nread,npuse,nouse,twind,gstime,sis,& mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i),read_rec(i),dval_use) @@ -1782,7 +1781,7 @@ subroutine read_obs(ndata,mype) ! Process ssmi data else if (obstype == 'ssmi' ) then - call read_ssmi(mype,val_dat,ithin,rmesh,platid,gstime,& + call read_ssmi(mype,val_dat,ithin,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis,& mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i),read_rec(i),dval_use) @@ -1791,7 +1790,7 @@ subroutine read_obs(ndata,mype) ! Process amsre data else if ( obstype == 'amsre_low' .or. obstype == 'amsre_mid' .or. & obstype == 'amsre_hig' ) then - call read_amsre(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,& + call read_amsre(mype,val_dat,ithin,isfcalc,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis,& mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i),read_rec(i),dval_use) @@ -1801,7 +1800,7 @@ subroutine read_obs(ndata,mype) else if (obstype == 'ssmis' .or. & obstype == 'ssmis_las' .or. obstype == 'ssmis_uas' .or. & obstype == 'ssmis_img' .or. obstype == 'ssmis_env' ) then - call read_ssmis(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,& + call read_ssmis(mype,val_dat,ithin,isfcalc,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis,& mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i),read_rec(i),dval_use) @@ -1809,7 +1808,7 @@ subroutine read_obs(ndata,mype) ! Process AMSR2 data else if(obstype == 'amsr2')then - call read_amsr2(mype,val_dat,ithin,rmesh,platid,gstime,& + call read_amsr2(mype,val_dat,ithin,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis,& mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i)) @@ -1817,7 +1816,7 @@ subroutine read_obs(ndata,mype) ! Process GOES IMAGER RADIANCE data else if(obstype == 'goes_img') then - call read_goesimg(mype,val_dat,ithin,rmesh,platid,gstime,& + call read_goesimg(mype,val_dat,ithin,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis, & mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i),read_rec(i),dval_use) @@ -1825,7 +1824,7 @@ subroutine read_obs(ndata,mype) ! Process GMI data else if (obstype == 'gmi') then - call read_gmi(mype,val_dat,ithin,rmesh,platid,gstime,& + call read_gmi(mype,val_dat,ithin,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis,& mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i),dval_use) @@ -1833,14 +1832,14 @@ subroutine read_obs(ndata,mype) ! Process Meteosat SEVIRI RADIANCE data else if(obstype == 'seviri') then - call read_seviri(mype,val_dat,ithin,rmesh,platid,gstime,& + call read_seviri(mype,val_dat,ithin,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis, & mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i),read_rec(i),dval_use) string='READ_SEVIRI' ! Process GOES-R ABI RADIANCE data else if(obstype == 'abi') then - call read_abi(mype,val_dat,ithin,rmesh,platid,gstime,& + call read_abi(mype,val_dat,ithin,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis, & mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i),read_rec(i),dval_use) @@ -1848,7 +1847,7 @@ subroutine read_obs(ndata,mype) ! Process Himawari-8 AHI RADIANCE data else if(obstype == 'ahi') then - call read_ahi(mype,val_dat,ithin,rmesh,platid,gstime,& + call read_ahi(mype,val_dat,ithin,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis, & mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i),read_rec(i),dval_use) @@ -1857,7 +1856,7 @@ subroutine read_obs(ndata,mype) ! Process NAVY AVHRR RADIANCE data else if(obstype == 'avhrr_navy') then - call read_avhrr_navy(mype,val_dat,ithin,rmesh,platid,gstime,& + call read_avhrr_navy(mype,val_dat,ithin,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis, & mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i),read_rec(i),dval_use) @@ -1865,14 +1864,14 @@ subroutine read_obs(ndata,mype) ! Process NESDIS AVHRR RADIANCE data else if(obstype == 'avhrr') then - call read_avhrr(mype,val_dat,ithin,rmesh,platid,gstime,& + call read_avhrr(mype,val_dat,ithin,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis, & mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i),read_rec(i),dval_use) string='READ_AVHRR' ! Process SST VIIRS RADIANCE data else if(obstype == 'viirs-m') then - call read_sst_viirs(mype,val_dat,ithin,rmesh,platid,gstime,& + call read_sst_viirs(mype,val_dat,ithin,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis, & mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i),read_rec(i),dval_use) @@ -1885,7 +1884,7 @@ subroutine read_obs(ndata,mype) if (is_extOzone(infile,obstype,dplat(i))) then call extOzone_read(infile,obstype,dplat(i),dsis(i), & - iread,ipuse,iouse, platid,gstime,lunout,twind,ithin,rmesh, & + iread,ipuse,iouse, dplat(i),gstime,lunout,twind,ithin,rmesh, & nobs_sub1(:,i)) string='extOzone_read' @@ -1895,7 +1894,7 @@ subroutine read_obs(ndata,mype) else call read_ozone(nread,npuse,nouse,& - platid,infile,gstime,lunout,obstype,twind,sis,ithin,rmesh, & + dplat(i),infile,gstime,lunout,obstype,twind,sis,ithin,rmesh, & nobs_sub1(1,i)) string='READ_OZONE' endif ozone_obstype_select @@ -1921,7 +1920,7 @@ subroutine read_obs(ndata,mype) ! Process aerosol data else if (ditype(i) == 'aero' )then call read_aerosol(nread,npuse,nouse,& - platid,infile,gstime,lunout,obstype,twind,sis,ithin,rmesh, & + dplat(i),infile,gstime,lunout,obstype,twind,sis,ithin,rmesh, & mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i)) string='READ_AEROSOL' @@ -1947,7 +1946,7 @@ subroutine read_obs(ndata,mype) call warn('read_obs',' infile =',trim(infile)) call warn('read_obs',' ditype(i) =',trim(ditype(i))) call warn('read_obs',' obstype =',trim(obstype)) - call warn('read_obs',' platid =',trim(platid)) + call warn('read_obs',' dplat =',trim(dplat(i))) call warn('read_obs',' string =',trim(string)) endif diff --git a/src/gsi/read_viirs.f90 b/src/gsi/read_viirs.f90 index ab80642f29..3ea50352ec 100644 --- a/src/gsi/read_viirs.f90 +++ b/src/gsi/read_viirs.f90 @@ -475,7 +475,7 @@ subroutine read_sst_viirs(mype,val_viirs,ithin,rmesh,jsatid,& nele,itxmax,nread,ndata_mesh,data_mesh,score_crit,nrec) if ( nread > 0 ) then - write(*,'(a,a10,I3,F6.1,3I10)') 'read_viirs,jsatid,imesh,amesh,itxmax,nread,ndata_mesh :',jsatid,imesh,amesh(imesh),itxmax,nread,ndata_mesh + write(*,'(a,a11,I3,F6.1,3I10)') 'read_viirs,jsatid,imesh,amesh,itxmax,nread,ndata_mesh :',jsatid,imesh,amesh(imesh),itxmax,nread,ndata_mesh endif ! ! get data_all by combining data from all thinning box sizes diff --git a/src/gsi/statsrad.f90 b/src/gsi/statsrad.f90 index d42a53f7d6..dfedcc92ab 100644 --- a/src/gsi/statsrad.f90 +++ b/src/gsi/statsrad.f90 @@ -100,7 +100,7 @@ subroutine statsrad(aivals,stats,ndata) write(iout_rad,4002) 'qcpenalty','qc1','qc2','qc3','qc4','qc5','qc6','qc7' write(iout_rad,4003) aivals(39,is),(nint(aivals(i,is)),i=8,14) 4000 format(2x,a3,7x,a4,14x,a7,1x,8(a7,1x)) -4001 format(1x,a10,1x,a10,f16.8,8(i7,1x)) +4001 format(1x,a11,1x,a10,f16.8,8(i7,1x)) 4002 format(28x,a9,1x,8(a7,1x)) 4003 format(22x,f16.8,8(i7,1x)) @@ -161,11 +161,11 @@ subroutine statsrad(aivals,stats,ndata) 2011 format(8x,f16.8,8(i7,1x)) 2012 format(12x,A7,5x,8(a7,1x)) 2999 format(' Illegal satellite type ') -1102 format(1x,i4,i5,1x,a16,2i7,1x,f10.3,1x,6(f11.7,1x)) +1102 format(1x,i4,i6,1x,a20,2i7,1x,f10.3,1x,6(f11.7,1x)) 1109 format(t5,'it',t13,'satellite',t23,'instrument',t40, & '# read',t53,'# keep',t65,'# assim',& t75,'penalty',t88,'qcpnlty',t104,'cpen',t115,'qccpen') -1115 format('o-g',1x,i2.2,1x,'rad',2x,2A10,2x,3(i11,2x),4(g12.5,1x)) +1115 format('o-g',1x,i2.2,1x,'rad',2x,2A12,2x,3(i11,2x),4(g12.5,1x)) ! Close output unit close(iout_rad)