From d72e074a48d6b080fcf724cd4425b3429dc469ee Mon Sep 17 00:00:00 2001 From: wx20jjung Date: Wed, 12 Jun 2024 12:19:19 +0000 Subject: [PATCH 01/16] Removed an unnecessary line --- src/gsi/read_cris.f90 | 1 - 1 file changed, 1 deletion(-) diff --git a/src/gsi/read_cris.f90 b/src/gsi/read_cris.f90 index d5668f8864..7e9195fe41 100644 --- a/src/gsi/read_cris.f90 +++ b/src/gsi/read_cris.f90 @@ -886,7 +886,6 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& 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) From 3c4147128709d519db4328dfdfea99fdc0afa067 Mon Sep 17 00:00:00 2001 From: wx20jjung Date: Wed, 12 Jun 2024 12:32:02 +0000 Subject: [PATCH 02/16] Initial coding for IASI-NG. Reads BUFR files generated by NESDIS. --- src/gsi/crtm_interface.f90 | 14 + src/gsi/gsi_files.cmake | 1 + src/gsi/gsimod.F90 | 11 +- src/gsi/obsmod.F90 | 2 +- src/gsi/qcmod.f90 | 42 +- src/gsi/read_iasing.f90 | 968 +++++++++++++++++++++++++++++++++++++ src/gsi/read_obs.F90 | 21 +- src/gsi/setuprad.f90 | 13 +- 8 files changed, 1043 insertions(+), 29 deletions(-) create mode 100644 src/gsi/read_iasing.f90 diff --git a/src/gsi/crtm_interface.f90 b/src/gsi/crtm_interface.f90 index 4bb1191001..8a7fb84507 100644 --- a/src/gsi/crtm_interface.f90 +++ b/src/gsi/crtm_interface.f90 @@ -662,6 +662,20 @@ subroutine init_crtm(init_pass,mype_diaghdr,mype,nchanl,nreal,isis,obstype,radmo error_status = crtm_channelinfo_subset(channelinfo(1), & channel_subset = nuchan(subset_start:subset_end)) +else if (channelinfo(1)%sensor_id(1:7) == 'iasi-ng' .AND. isis(1:7) == 'iasi-ng') then + sensorindex = 1 + subset_start = 0 + subset_end = 0 + do k=1, jpch_rad + if (isis == nusis(k)) then + if (subset_start == 0) subset_start = k + subset_end = k + endif + end do + + error_status = crtm_channelinfo_subset(channelinfo(1), & + channel_subset = nuchan(subset_start:subset_end)) + else if (channelinfo(1)%sensor_id(1:4) == 'iasi' .AND. isis(1:4) == 'iasi') then sensorindex = 1 subset_start = 0 diff --git a/src/gsi/gsi_files.cmake b/src/gsi/gsi_files.cmake index 95d885e2ee..7a08d9e847 100644 --- a/src/gsi/gsi_files.cmake +++ b/src/gsi/gsi_files.cmake @@ -493,6 +493,7 @@ read_goesndr.f90 read_gps.f90 read_guess.F90 read_iasi.f90 +read_iasing.f90 read_l2bufr_mod.f90 read_lag.f90 read_lidar.f90 diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index d7f5667252..d11616b55d 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -94,7 +94,7 @@ module gsimod init_qcvars,vadfile,noiqc,c_varqc,gps_jacqc,qc_noirjaco3,qc_noirjaco3_pole,& buddycheck_t,buddydiag_save,njqc,vqc,nvqc,hub_norm,vadwnd_l2rw_qc, & pvis,pcldch,scale_cv,estvisoe,estcldchoe,vis_thres,cldch_thres,cao_check, & - cris_cads, iasi_cads, airs_cads + cris_cads, iasi_cads, iasing_cads, airs_cads use qcmod, only: troflg,lat_c,nrand use cads, only: M__Sensor,N__Num_Bands,N__GradChkInterval,N__Band_Size,N__Bands,N__Window_Width, & N__Window_Bounds,R__BT_Threshold,R__Grad_Threshold,R__Window_Grad_Threshold, L__Do_Quick_Exit, & @@ -1064,9 +1064,10 @@ module gsimod ! ! Flags to use the new IR cloud detection routine. Flag must be set to true to use the new routine. The default ! (no flag or .false.) will use the default. -! airs_cads: use the clod and aerosool detection software for the AIRS instrument -! cris_cads: use the cloud and aerosol detection software for CrIS instruments -! iasi_cads: use the cloud and aerosol detection software for IASI instruments +! airs_cads : use the clod and aerosool detection software for the AIRS instrument +! cris_cads : use the cloud and aerosol detection software for CrIS instruments +! iasi_cads : use the cloud and aerosol detection software for IASI instruments +! iasing_cads: use the cloud and aerosol detection software for IASI instruments ! namelist/obsqc/dfact,dfact1,erradar_inflate,tdrerr_inflate,oberrflg,& @@ -1078,7 +1079,7 @@ module gsimod q_doe_a_136,q_doe_a_137,q_doe_b_136,q_doe_b_137, & t_doe_a_136,t_doe_a_137,t_doe_b_136,t_doe_b_137, & uv_doe_a_236,uv_doe_a_237,uv_doe_a_213,uv_doe_b_236,uv_doe_b_237,uv_doe_b_213, & - vad_near_analtime,airs_cads,cris_cads,iasi_cads + vad_near_analtime,airs_cads,cris_cads,iasi_cads, iasing_cads ! OBS_INPUT (controls input data): ! dmesh(max(dthin))- thinning mesh for each group diff --git a/src/gsi/obsmod.F90 b/src/gsi/obsmod.F90 index 1c45c62bc8..95e5046d9f 100644 --- a/src/gsi/obsmod.F90 +++ b/src/gsi/obsmod.F90 @@ -1182,7 +1182,7 @@ subroutine init_obsmod_vars(nhr_assim,mype) endif ! for cris, iasi, atms, regional analysis may want shorter time window if (index(dtype(ii),'cris') /= 0 .or. index(dtype(ii),'atms') /= 0 .or. & - index(dtype(ii),'iasi') /= 0 ) then + index(dtype(ii),'iasi') /= 0 .or. index(dtype(ii),'iasi-ng') /= 0) then if(time_window(ii)>time_window_rad) then time_window(ii) = time_window_rad if (mype==0) write(6,*) 'INIT_OBSMOD_VARS: reset time window for ',dtype(ii),& diff --git a/src/gsi/qcmod.f90 b/src/gsi/qcmod.f90 index 9804965573..6425d5268e 100644 --- a/src/gsi/qcmod.f90 +++ b/src/gsi/qcmod.f90 @@ -118,6 +118,7 @@ module qcmod ! def airs_cads - if true, use the cloud and aerosol detection routine for Aqua/AIRS instrument ! def cris_cads - if true, use the cloud and aerosol detection routine for CrIS instruments ! def iasi_cads - if true, use the cloud and aerosol detection routine for IASI instruments +! def iasing_cads - if true, use the cloud and aerosol detection routine for IASI-NG instruments ! ! following used for nonlinear qc: ! @@ -204,7 +205,7 @@ module qcmod public :: troflg public :: lat_c public :: nrand - public :: airs_cads, cris_cads, iasi_cads + public :: airs_cads, cris_cads, iasi_cads, iasing_cads logical nlnqc_iter,njqc,vqc,nvqc,hub_norm logical noiqc @@ -220,7 +221,7 @@ module qcmod logical vadwnd_l2rw_qc logical troflg logical cao_check - logical airs_cads, cris_cads, iasi_cads + logical airs_cads, cris_cads, iasi_cads, iasing_cads character(10):: vadfile integer(i_kind) npres_print @@ -461,9 +462,10 @@ subroutine init_qcvars lat_c=21.0_r_kind nrand=13 - airs_cads = .false. - cris_cads = .false. - iasi_cads = .false. + airs_cads = .false. + cris_cads = .false. + iasi_cads = .false. + iasing_cads = .false. return end subroutine init_qcvars @@ -598,8 +600,8 @@ subroutine setup_tzr_qc(obstype) tzchk = 0.85_r_kind elseif ( obstype == 'hirs2' .or. obstype == 'hirs3' .or. obstype == 'hirs4' .or. & obstype == 'sndr' .or. obstype == 'sndrd1' .or. obstype == 'sndrd2'.or. & - obstype == 'sndrd3' .or. obstype == 'sndrd4' .or. & - obstype == 'goes_img' .or. obstype == 'ahi' .or. obstype == 'airs' .or. obstype == 'iasi' .or. & + obstype == 'sndrd3' .or. obstype == 'sndrd4' .or. obstype == 'goes_img' .or. & + obstype == 'ahi' .or. obstype == 'airs' .or. obstype == 'iasi' .or. obstype == 'iasi-ng' .or.& obstype == 'cris' .or. obstype == 'cris-fsr' .or. obstype == 'seviri' .or. obstype == 'abi') then tzchk = 0.85_r_kind endif @@ -2076,7 +2078,7 @@ subroutine qc_saphir(nchanl,sfchgt,luse,sea, & end subroutine qc_saphir subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr,airs, & - cris,iasi,hirs,zsges,cenlat,frac_sea,pangs,trop5,zasat,tzbgr,tsavg5,tbc,tb_obs,tbcnob,tnoise, & + cris,iasi,iasing,hirs,zsges,cenlat,frac_sea,pangs,trop5,zasat,tzbgr,tsavg5,tbc,tb_obs,tbcnob,tnoise, & wavenumber,ptau5,prsltmp,tvp,temp,wmix,chan_level,emissivity_k,ts,tsim, & id_qc,aivals,errf,varinv,varinv_use,cld,cldp,kmax,zero_irjaco3_pole,cluster_fraction, & cluster_bt, chan_stdev, model_bt) @@ -2084,12 +2086,12 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr,airs !$$$ subprogram documentation block ! . . . -! subprogram: qc_irsnd QC for ir sounder data(hirs,goessndr,airs,iasi,cris) +! subprogram: qc_irsnd QC for ir sounder data(hirs,goessndr,airs,iasi,iasing,cris) ! ! prgmmr: derber org: np23 date: 2010-08-20 ! ! abstract: set quality control criteria for ir sounder data (hirs, -! goessndr, airs, iasi, cris) +! goessndr, airs, iasi, iasing, cris) ! ! program history log: ! 2010-08-10 derber transfered from setuprad @@ -2164,7 +2166,7 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr,airs ! Declare passed variables - logical, intent(in ) :: sea,land,ice,snow,luse,goessndr,airs,cris,hirs,iasi + logical, intent(in ) :: sea,land,ice,snow,luse,goessndr,airs,cris,hirs,iasi,iasing logical, intent(inout) :: zero_irjaco3_pole integer(i_kind), intent(in ) :: nsig,nchanl,ndat,is integer(i_kind),dimension(nchanl), intent(in ) :: ich @@ -2326,6 +2328,21 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr,airs tropopause_height, boundary_layer_pres, tb_bc, tsim, chan_level, imager_chans, cluster_fraction, & cluster_bt, chan_stdev, model_bt, i_flag_cloud, cldp ) + elseif ( iasing .and. iasing_cads ) then + I_Sensor_ID = 59 + chan_array = nuchan(ich) ! channel numbers + tb_bc = tbc + tsim ! observation BT with bias correction + boundary_layer_pres = nint(0.8_r_kind*prsltmp(1)) ! boundary layer set to be 80% of surface pressure + tropopause_height = nint(trop5) + imager_chans = (/18,19/) ! imager channel numbers (from satinfo) + isurface_chan = 2447 ! surface channel + ichan_10_micron = 2447 ! ~10.7 micron channel for low level cloud test + ichan_12_micron = 1560 ! ~12.0 micron channel for low level cloud test + + call cloud_aerosol_detection( I_Sensor_ID, nchanl, chan_array, & + tropopause_height, boundary_layer_pres, tb_bc, tsim, chan_level, imager_chans, cluster_fraction, & + cluster_bt, chan_stdev, model_bt, i_flag_cloud, cldp ) + elseif ( airs .and. airs_cads ) then I_Sensor_ID = 11 chan_array = nuchan(ich) ! channel numbers @@ -2348,7 +2365,8 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr,airs ! compute cloud stats ! If using CADS - if ((cris .and. cris_cads) .or. (iasi .and. iasi_cads) .or. (airs .and. airs_cads)) then + if ((cris .and. cris_cads) .or. (iasi .and. iasi_cads) .or. (airs .and. airs_cads) .or. & + iasing .and. iasing_cads ) then ! Reject channels affected by clouds do i=1, nchanl diff --git a/src/gsi/read_iasing.f90 b/src/gsi/read_iasing.f90 new file mode 100644 index 0000000000..d0ebb9d9c1 --- /dev/null +++ b/src/gsi/read_iasing.f90 @@ -0,0 +1,968 @@ +subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& + infile,lunout,obstype,nread,ndata,nodata,twind,sis,& + mype_root,mype_sub,npe_sub,mpi_comm_sub,nobs, & + nrec_start,nrec_start_ears,nrec_start_db,dval_use) +!$$$ subprogram documentation block +! . . . . +! subprogram: read_iasing read bufr format iasing data +! prgmmr : jung org: date: 2024-05-30 +! +! abstract: This routine reads BUFR format radiance +! files. Optionally, the data are thinned to +! a specified resolution using simple quality control checks. +! +! When running the gsi in regional mode, the code only +! retains those observations that fall within the regional +! domain +! +! program history log: +! 2024-06-12 jung - transposed read_iasi subroutine to read_iasing +! +! +! +! input argument list: +! mype - mpi task id +! val_iasing - weighting factor applied to super obs +! ithin - flag to thin data +! isfcalc - when set to one, calculate surface characteristics using +! method that accounts for the size/shape of the fov. +! when not one, calculate surface characteristics using +! bilinear interpolation. +! rmesh - thinning mesh size (km) +! jsatid - satellite id +! gstime - analysis time in minutes from reference date +! infile - unit from which to read BUFR data +! lunout - unit to which to write data for further processing +! obstype - observation type to process +! twind - input group time window (hours) +! sis - sensor/instrument/satellite indicator +! mype_root - "root" task for sub-communicator +! mype_sub - mpi task id within sub-communicator +! npe_sub - number of data read tasks +! mpi_comm_sub - sub-communicator for data read +! nrec_start - first subset with useful information +! nrec_start_ears - first ears subset with useful information +! nrec_start_db - first db subset with useful information +! +! output argument list: +! nread - number of BUFR IASI-NG observations read +! ndata - number of BUFR IASI-NG profiles retained for further processing +! nodata - number of BUFR IASI-NG observations retained for further processing +! nobs - array of observations on each subdomain for each processor +! +! attributes: +! language: f90 +! machine: ibm RS/6000 SP +! +!$$$ +! Use modules + use kinds, only: r_kind,r_double,i_kind + use satthin, only: super_val,itxmax,makegrids,map2tgrid,destroygrids, & + finalcheck,checkob,score_crit + use satthin, only: radthin_time_info,tdiff2crit + use obsmod, only: time_window_max + use radinfo, only:iuse_rad,nuchan,nusis,jpch_rad,crtm_coeffs_path,use_edges, & + radedge1,radedge2,radstart,radstep + use crtm_module, only: success, & + crtm_kind => fp + use crtm_planck_functions, only: crtm_planck_temperature + use crtm_spccoeff, only: sc,crtm_spccoeff_load,crtm_spccoeff_destroy + use gridmod, only: diagnostic_reg,regional,nlat,nlon,& + tll2xy,txy2ll,rlats,rlons + use constants, only: zero,deg2rad,rad2deg,r60inv,one,ten,r100 + use gsi_4dvar, only: l4dvar,l4densvar,iwinbgn,winlen + use calc_fov_crosstrk, only: instrument_init, fov_check, fov_cleanup + use deter_sfc_mod, only: deter_sfc,deter_sfc_fov + use obsmod, only: bmiss + use gsi_nstcouplermod, only:nst_gsi,nstinfo + use gsi_nstcouplermod, only: gsi_nstcoupler_skindepth, gsi_nstcoupler_deter + use mpimod, only: npe + use gsi_io, only: verbose + use qcmod, only: iasing_cads +! use radiance_mod, only: rad_obs_type + + implicit none + +! BUFR format for IASISPOT +! Input variables + integer(i_kind) ,intent(in ) :: mype,nrec_start,nrec_start_ears,nrec_start_db + integer(i_kind) ,intent(in ) :: ithin + integer(i_kind) ,intent(inout) :: isfcalc + integer(i_kind) ,intent(in ) :: lunout + integer(i_kind) ,intent(in ) :: mype_root + 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=20),intent(in ) :: sis + real(r_kind) ,intent(in ) :: twind + real(r_kind) ,intent(inout) :: val_iasing + real(r_kind) ,intent(in ) :: gstime + real(r_kind) ,intent(in ) :: rmesh + logical ,intent(in ) :: dval_use + +! Output variables + integer(i_kind) ,intent(inout) :: nread + integer(i_kind),dimension(npe) ,intent(inout) :: nobs + integer(i_kind) ,intent( out) :: ndata,nodata + + +! BUFR file sequencial number +! character(len=512) :: table_file + integer(i_kind) :: lnbufr = 10 + +! Variables for BUFR IO + real(r_double) :: crchn_reps + real(r_double),dimension(5) :: linele + real(r_double),dimension(13) :: allspot + real(r_double),allocatable,dimension(:,:) :: allchan + real(r_double),dimension(3,4):: cscale + real(r_double),dimension(7):: cloud_frac + integer(i_kind) :: bufr_size + + real(r_kind) :: step, start,step_adjust + character(len=8) :: subset + character(len=4) :: senname + character(len=80) :: allspotlist + character(len=40) :: infile2 + integer(i_kind) :: iret,ireadsb,ireadmg,irec,next, nrec_startx + integer(i_kind),allocatable,dimension(:) :: nrec + + +! Work variables for time + integer(i_kind) :: idate + integer(i_kind) :: idate5(5) + real(r_kind) :: sstime, tdiff, t4dv + integer(i_kind) :: nmind + + +! Other work variables + real(r_kind) :: piece + real(r_kind) :: rsat, dlon, dlat + real(r_kind) :: dlon_earth,dlat_earth,dlon_earth_deg,dlat_earth_deg + real(r_kind) :: lza, lzaest,sat_height_ratio + real(r_kind) :: pred, crit1, dist1 + real(r_kind) :: sat_zenang + real(crtm_kind) :: radiance + real(r_kind) :: tsavg,vty,vfr,sty,stp,sm,sn,zz,ff10,sfcr + real(r_kind) :: zob,tref,dtw,dtc,tz_tr + real(r_kind),dimension(0:4) :: rlndsea + real(r_kind),dimension(0:3) :: sfcpct + real(r_kind),dimension(0:3) :: ts + real(r_kind),dimension(10) :: sscale + real(crtm_kind),allocatable,dimension(:) :: temperature + real(r_kind),allocatable,dimension(:) :: scalef + real(r_kind),allocatable,dimension(:,:):: data_all + real(r_kind) cdist,disterr,disterrmax,dlon00,dlat00 + + logical :: outside,iuse,assim,valid + logical :: iasing,quiet,cloud_info + + integer(i_kind) :: ifov, instr, iscn, ioff, sensorindex_iasing + integer(i_kind) :: i, j, l, iskip, ifovn, bad_line, ksatid, kidsat, llll + integer(i_kind) :: nreal, isflg + integer(i_kind) :: itx, k, nele, itt, n + integer(i_kind):: iexponent,maxinfo, bufr_nchan, dval_info + integer(i_kind):: idomsfc(1) + integer(i_kind):: ntest + integer(i_kind):: error_status, irecx,ierr + integer(i_kind):: radedge_min, radedge_max + integer(i_kind) :: subset_start, subset_end, satinfo_nchan, sc_chan, bufr_chan + integer(i_kind) :: sfc_channel_index + integer(i_kind),allocatable, dimension(:) :: channel_number, sc_index, bufr_index + integer(i_kind),allocatable, dimension(:) :: bufr_chan_test + character(len=20),allocatable, dimension(:):: sensorlist + +! Imager cluster information for CADS + integer(i_kind) :: sensorindex_imager, cads_info + integer(i_kind),dimension(7) :: imager_cluster_index + logical :: imager_coeff + logical,dimension(7) :: imager_cluster_flag + character(len=80) :: spc_filename + real(r_kind),dimension(123,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 + +! Set standard parameters + character(8),parameter:: fov_flag="crosstrk" + integer(i_kind),parameter:: sfc_channel=1271 + integer(i_kind),parameter:: ichan=-999 ! fov-based surface code is not channel specific for iasing + real(r_kind),parameter:: expansion=one ! exansion factor for fov-based surface code. + ! use one for ir sensors. + real(r_kind),parameter:: R90 = 90._r_kind + real(r_kind),parameter:: R360 = 360._r_kind + real(r_kind),parameter:: tbmin = 50._r_kind + real(r_kind),parameter:: tbmax = 550._r_kind + real(r_kind),parameter:: earth_radius = 6371000._r_kind + integer(i_kind),parameter :: ilon = 3 + integer(i_kind),parameter :: ilat = 4 + real(r_kind) :: ptime,timeinflat,crit0 + integer(i_kind) :: ithin_time,n_tbin,it_mesh,jstart + logical print_verbose + + print_verbose=.false. + if(verbose)print_verbose=.true. + +! Initialize variables + maxinfo = 31 + disterrmax=zero + ntest=0 + dval_info = 0 + if(dval_use) dval_info = 2 + cads_info = 0 + if(iasing_cads) cads_info = 23 + nreal = maxinfo + cads_info + dval_info + nstinfo + + ndata = 0 + nodata = 0 + iasing = obstype == 'iasi-ng' + + bad_line=-1 + + if (nst_gsi > 0 ) then + call gsi_nstcoupler_skindepth(obstype, zob) ! get penetration depth (zob) for the obstype + endif + + if(jsatid == 'metop-sg-a1')kidsat=24 + +! write(6,*)'READ_IASI-NG: mype, mype_root,mype_sub, npe_sub,mpi_comm_sub', & +! mype, mype_root,mype_sub,mpi_comm_sub + + radedge_min = 0 + radedge_max = 1000 + +! Find the iasi-ng offset in the jpch_rad list. This is for the iuse flag +! and count the number of cahnnels in the satinfo file + ioff=jpch_rad + subset_start = 0 + subset_end = 0 + assim = .false. + do i=1,jpch_rad + if (trim(nusis(i))==trim(sis)) then + ioff = min(ioff,i) ! iasi-ng offset + if (subset_start == 0) then + step = radstep(i) + start = radstart(i) + if (radedge1(i)/=-1 .and. radedge2(i)/=-1) then + radedge_min=radedge1(i) + radedge_max=radedge2(i) + end if + subset_start = i + endif + if (iuse_rad(i) > 0) assim = .true. ! Are any of the IASI-NG channels being used? + subset_end = i + endif + end do + satinfo_nchan = subset_end - subset_start + 1 + allocate(channel_number(satinfo_nchan)) + allocate(sc_index(satinfo_nchan)) + allocate(bufr_index(satinfo_nchan)) + ioff = ioff - 1 + + step_adjust = 0.625_r_kind +! If all channels of a given sensor are set to monitor or not +! assimilate mode (iuse_rad<1), reset relative weight to zero. +! We do not want such observations affecting the relative +! weighting between observations within a given thinning group. + if (.not.assim) val_iasing=zero + + if (mype_sub==mype_root)write(6,*)'READ_IASI-NG: iasi-ng offset ',ioff + + senname = 'IASI-NG' + + allspotlist= & + 'SAID YEAR MNTH DAYS HOUR MINU SECO CLATH CLONH SAZA BEARAZ SOZA SOLAZI' + +! load spectral coefficient structure + quiet=.not. verbose + + imager_coeff = .false. + spc_filename =trim(crtm_coeffs_path)//'metimage_'//trim(jsatid)//'.SpcCoeff.bin' + inquire(file=trim(spc_filename), exist=imager_coeff) + if ( imager_coeff ) then + allocate( sensorlist(2)) + sensorlist(1) = sis + sensorlist(2) = 'metimage_'//trim(jsatid) + else + allocate( sensorlist(1)) + sensorlist(1) = sis + endif + + if( crtm_coeffs_path /= "" ) then + if(mype_sub==mype_root .and. print_verbose) write(6,*)'READ_IASI-NG: crtm_spccoeff_load() on path "'//trim(crtm_coeffs_path)//'"' + error_status = crtm_spccoeff_load(sensorlist,& + File_Path = crtm_coeffs_path,quiet=quiet ) + else + error_status = crtm_spccoeff_load(sensorlist,quiet=quiet) + endif + + if (error_status /= success) then + write(6,*)'READ_IASI-NG: ***ERROR*** crtm_spccoeff_load error_status=',error_status,& + ' TERMINATE PROGRAM EXECUTION' + call stop2(71) + endif + +! find IASI-NG sensorindex + sensorindex_iasing = 0 + if ( sc(1)%sensor_id(1:7) == 'iasi-ng' ) then + sensorindex_iasing = 1 + else + write(6,*)'READ_IASI-NG: ***ERROR*** sensorindex_iasi-ng not set NO IASI-NG DATA USED' + write(6,*)'READ_IASI-NG: We are looking for ', sc(1)%sensor_id, ' TERMINATE PROGRAM EXECUTION' + call stop2(71) + end if + +! find imager sensorindex + sensorindex_imager = 0 + if ( iasing_cads .and. imager_coeff ) then + if ( sc(2)%sensor_id(1:7) == 'metimag' ) then + sensorindex_imager = 2 + imager_coeff = .true. + else + write(6,*)'READ_IASI-NG: ***ERROR*** sensorindex_imager is not set' + write(6,*)'READ_IASI-NG: We are looking for ', sc(2)%sensor_id + imager_coeff = .false. + end if + else + imager_coeff = .false. + end if + +! Find the channels being used (from satinfo file) in the spectral coef. structure. + do i=subset_start,subset_end + channel_number(i -subset_start +1) = nuchan(i) + end do + sc_index(:) = 0 + satinfo_chan: do i=1,satinfo_nchan + spec_coef: do l=1,sc(1)%n_channels + if ( channel_number(i) == sc(sensorindex_iasing)%sensor_channel(l) ) then + sc_index(i) = l + exit spec_coef + endif + end do spec_coef + end do satinfo_chan + +! Calculate parameters needed for FOV-based surface calculation. + if (isfcalc==1)then + instr=18 + call instrument_init(instr, jsatid, expansion, valid) + if (.not. valid) then + if (assim) then + write(6,*)'READ_IASI-NG: ***ERROR*** IN SETUP OF FOV-SFC CODE. STOP' + call stop2(71) + else + call fov_cleanup + isfcalc = 0 + write(6,*)'READ_IASI-NG: ***ERROR*** IN SETUP OF FOV-SFC CODE' + endif + endif + endif + + if (isfcalc==1)then + rlndsea = zero + else + rlndsea(0) = zero + rlndsea(1) = 10._r_kind + rlndsea(2) = 15._r_kind + rlndsea(3) = 10._r_kind + rlndsea(4) = 30._r_kind + endif + + call radthin_time_info(obstype, jsatid, sis, ptime, ithin_time) + if( ptime > 0.0_r_kind) then + n_tbin=nint(2*time_window_max/ptime) + else + n_tbin=1 + endif +! Make thinning grids + call makegrids(rmesh,ithin,n_tbin=n_tbin) + +! Allocate arrays to hold data +! The number of channels is obtained from the satinfo file being used. + nele=nreal+satinfo_nchan + allocate(data_all(nele,itxmax),nrec(itxmax)) + allocate(temperature(1)) ! dependent on # of channels in the bufr file + allocate(allchan(2,1)) ! actual values set after ireadsb + allocate(bufr_chan_test(1))! actual values set after ireadsb + allocate(scalef(1)) + +! Big loop to read data file + next=0 + irec=0 + nrec=999999 +! Big loop over standard data feed and possible rars/db data +! llll=1 is normal feed, llll=2 RARS/EARS data, llll=3 DB/UW data) + ears_db_loop: do llll= 1, 3 + + if(llll == 1)then + nrec_startx=nrec_start + infile2=trim(infile) ! Set bufr subset names based on type of data to read + elseif(llll == 2) then + nrec_startx=nrec_start_ears + infile2=trim(infile)//'ears' ! Set bufr subset names based on type of data to read + elseif(llll == 3) then + nrec_startx=nrec_start_db + infile2=trim(infile)//'_db' ! Set bufr subset names based on type of data to read + end if + +! Open BUFR file + open(lnbufr,file=trim(infile2),form='unformatted',status='old',iostat=ierr) + + if(ierr /= 0) cycle ears_db_loop +! Open BUFR table + call openbf(lnbufr,'IN',lnbufr) + call datelen(10) + + irecx = 0 + read_subset: do while(ireadmg(lnbufr,subset,idate)>=0) + irecx = irecx + 1 + if(irecx < nrec_startx) cycle read_subset + irec = irec + 1 + next=next+1 + if(next == npe_sub)next=0 + if(next /= mype_sub)cycle read_subset + + read_loop: do while (ireadsb(lnbufr)==0) + +! Get the size of the channels and radiance (allchan) array + call ufbint(lnbufr,crchn_reps,1,1,iret,'(I1CRSQ)') + bufr_nchan = int(crchn_reps) + + bufr_size = size(temperature,1) + if ( bufr_size /= bufr_nchan ) then ! Re-allocation if number of channels has changed +! Allocate the arrays needed for the channel and radiance array + deallocate(temperature,allchan,bufr_chan_test,scalef) + allocate(temperature(bufr_nchan)) ! dependent on # of channels in the bufr file + allocate(allchan(2,bufr_nchan)) + allocate(bufr_chan_test(bufr_nchan)) + allocate(scalef(bufr_nchan)) + bufr_chan_test(:)=0 + endif ! allocation if + +! Read IASI-NG FOV information + call ufbint(lnbufr,linele,5,1,iret,'FOVN SLNM INGGQF SELV SAID') + +! Extract satellite id. If not the one we want, read next subset + ksatid=nint(linele(5)) + if(ksatid /= kidsat) cycle read_loop + + if ( linele(3) /= zero) cycle read_loop ! problem with profile (INGGQF) + +! zenith angle/scan spot mismatch, reject entire line + if ( bad_line == nint(linele(2))) then + cycle read_loop + else + bad_line = -1 + endif + + ifov = nint(linele(1)) ! field of view + +! IASI-NG fov ranges from 1 to 120. Current angle dependent bias +! correction has a maximum of 90 scan positions. Geometry +! of IASI-NG scan allows us to remap 1-120 to 1-60. Variable +! ifovn below contains the remapped IASI-NG fov. This value is +! passed on to and used in setuprad + ifovn = (ifov-1)/2 + 1 + +! Remove data on edges + if (.not. use_edges .and. & + (ifovn < radedge_min .OR. ifovn > radedge_max )) cycle read_loop + +! Check field of view (FOVN) and satellite zenith angle (SAZA) + iscn = nint(linele(2)) ! scan line + if( ifov <= 0 .or. ifov > 120) then + write(6,*)'READ_IASI-NG: ### ERROR IN READING ', senname, ' BUFR DATA:', & + ' STRANGE OBS INFO(FOVN,SLNM):', ifov, iscn + cycle read_loop + endif + + call ufbint(lnbufr,allspot,13,1,iret,allspotlist) + if(iret /= 1) cycle read_loop + + +! Check observing position + dlat_earth = allspot(8) ! latitude + dlon_earth = allspot(9) ! longitude + if( abs(dlat_earth) > R90 .or. abs(dlon_earth) > R360 .or. & + (abs(dlat_earth) == R90 .and. dlon_earth /= ZERO) )then + write(6,*)'READ_IASI-NG: ### ERROR IN READING ', senname, ' BUFR DATA:', & + ' STRANGE OBS POINT (LAT,LON):', dlat_earth, dlon_earth + cycle read_loop + endif + +! Retrieve observing position + if(dlon_earth >= R360)then + dlon_earth = dlon_earth - R360 + else if(dlon_earth < ZERO)then + dlon_earth = dlon_earth + R360 + endif + + dlat_earth_deg = dlat_earth + dlon_earth_deg = dlon_earth + dlat_earth = dlat_earth * deg2rad + dlon_earth = dlon_earth * deg2rad + +! If regional, map obs lat,lon to rotated grid. + if(regional)then + +! Convert to rotated coordinate. dlon centered on 180 (pi), +! so always positive for limited area + call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside) + if(diagnostic_reg) then + call txy2ll(dlon,dlat,dlon00,dlat00) + ntest=ntest+1 + cdist=sin(dlat_earth)*sin(dlat00)+cos(dlat_earth)*cos(dlat00)* & + (sin(dlon_earth)*sin(dlon00)+cos(dlon_earth)*cos(dlon00)) + cdist=max(-one,min(cdist,one)) + disterr=acos(cdist)*rad2deg + disterrmax=max(disterrmax,disterr) + end if + +! Check to see if in domain. outside=.true. if dlon_earth, +! dlat_earth outside domain, =.false. if inside + if(outside) cycle read_loop + +! Global case + else + dlat = dlat_earth + dlon = dlon_earth + call grdcrd1(dlat,rlats,nlat,1) + call grdcrd1(dlon,rlons,nlon,1) + endif + +! Check obs time + idate5(1) = nint(allspot(2)) ! year + idate5(2) = nint(allspot(3)) ! month + idate5(3) = nint(allspot(4)) ! day + idate5(4) = nint(allspot(5)) ! hour + idate5(5) = nint(allspot(6)) ! minute + + if( idate5(1) < 1900 .or. idate5(1) > 3000 .or. & + idate5(2) < 1 .or. idate5(2) > 12 .or. & + idate5(3) < 1 .or. idate5(3) > 31 .or. & + idate5(4) < 0 .or. idate5(4) > 24 .or. & + idate5(5) < 0 .or. idate5(5) > 60 )then + + write(6,*)'READ_IASI-NG: ### ERROR IN READING ', senname, ' BUFR DATA:', & + ' STRANGE OBS TIME (YMDHM):', idate5(1:5) + cycle read_loop + + endif + +! Retrieve obs time + call w3fs21(idate5,nmind) + t4dv = (real(nmind-iwinbgn,r_kind) + real(allspot(7),r_kind)*r60inv)*r60inv ! add in seconds + sstime = real(nmind,r_kind) + real(allspot(7),r_kind)*r60inv ! add in seconds + tdiff = (sstime - gstime)*r60inv + + if (l4dvar.or.l4densvar) then + if (t4dvwinlen) cycle read_loop + else + if (abs(tdiff)>twind) cycle read_loop + endif + +! Increment nread counter by satinfo_nchan + nread = nread + satinfo_nchan + + crit0 = 0.01_r_kind + if( llll > 1 ) crit0 = crit0 + r100 * real(llll,r_kind) + timeinflat=6.0_r_kind + call tdiff2crit(tdiff,ptime,ithin_time,timeinflat,crit0,crit1,it_mesh) + call map2tgrid(dlat_earth,dlon_earth,dist1,crit1,itx,ithin,itt,iuse,sis,it_mesh=it_mesh) + + if(.not. iuse)cycle read_loop + +! Observational info + sat_zenang = allspot(10) ! satellite zenith angle + +! Check satellite zenith angle (SAZA) + if(sat_zenang > 90._r_kind ) then + write(6,*)'READ_IASI-NG: ### ERROR IN READING ', senname, ' BUFR DATA:', & + ' STRANGE OBS INFO(FOVN,SLNM,SAZA,BEARAZ):', ifov, iscn, allspot(10),allspot(11) + cycle read_loop + endif + if ( ifov <= 60 ) sat_zenang = -sat_zenang + +! Compare IASI-NG satellite scan angle and zenith angle + piece = -step_adjust + if ( mod(ifovn,2) == 1) piece = step_adjust + lza = ((start + real((ifov-1)/4,r_kind)*step) + piece)*deg2rad + sat_height_ratio = (earth_radius + linele(4))/earth_radius + lzaest = asin(sat_height_ratio*sin(lza))*rad2deg + if (abs(sat_zenang - lzaest) > one) then + write(6,*)' READ_IASI-NG WARNING uncertainty in lza ', & + lza*rad2deg,sat_zenang,sis,ifov,start,step,allspot(11),allspot(12),allspot(13) + bad_line = iscn + cycle read_loop + endif + + +! "Score" observation. We use this information to identify "best" obs +! Locate the observation on the analysis grid. Get sst and land/sea/ice +! mask. +! isflg - surface flag +! 0 sea +! 1 land +! 2 sea ice +! 3 snow +! 4 mixed + +! When using FOV-based surface code, must screen out obs with bad fov numbers. + if (isfcalc == 1) then + call fov_check(ifov,instr,ichan,valid) + if (.not. valid) cycle read_loop + +! When isfcalc is set to one, calculate surface fields using size/shape of fov. +! Otherwise, use bilinear interpolation. + + call deter_sfc_fov(fov_flag,ifov,instr,ichan,real(allspot(11),r_kind),dlat_earth_deg, & + dlon_earth_deg,expansion,t4dv,isflg,idomsfc(1), & + sfcpct,vfr,sty,vty,stp,sm,ff10,sfcr,zz,sn,ts,tsavg) + else + call deter_sfc(dlat,dlon,dlat_earth,dlon_earth,t4dv,isflg,idomsfc(1),sfcpct, & + ts,tsavg,vty,vfr,sty,stp,sm,sn,zz,ff10,sfcr) + endif + +! Set common predictor parameters + crit1 = crit1 + rlndsea(isflg) + + call checkob(one,crit1,itx,iuse) + if(.not. iuse)cycle read_loop + +! Clear Amount (percent clear) +! Compute "score" for observation. All scores>=0.0. Lowest score is "best" + pred = r100 + cloud_info = .false. + call ufbrep(lnbufr,cloud_frac,1,7,iret,'ASCS') + if (iret == 7 .and. cloud_frac(1) <= r100 .and. cloud_frac(1) >= zero) then + pred = r100 - cloud_frac(1) + cloud_info = .true. + endif + + crit1 = crit1 + pred + call checkob(one,crit1,itx,iuse) + if(.not. iuse)cycle read_loop + + call ufbseq(lnbufr,cscale,3,4,iret,'IAS1CBSQ') + if(iret /= 4) then + write(6,*) 'READ_IASI-NG read scale error ',iret + cycle read_loop + end if + +! The scaling factors are as follows, cscale(1) is the start channel number, +! cscale(2) is the end channel number, +! cscale(3) is the exponent scaling factor +! In our case, there are 4 groups of cscale (dimension :: cscale(3,4)) +! The units are W/m2..... you need to convert to mW/m2.... (subtract 5 from cscale(3,:) + do i=1,4 ! convert exponent scale factor to int and change units + if(cscale(3,i) < bmiss) then + iexponent = -(nint(cscale(3,i)) - 5) + sscale(i)=ten**iexponent + else + sscale(i)=0.0_r_kind + endif + end do + +!JAJ NEED TO CHECK THIS ONB +! Read IASI-NG channel number(CHNM) and radiance (SCRA) + call ufbseq(lnbufr,allchan,2,bufr_nchan,iret,'I1CRSQ') + jstart=1 + scalef=one + do i=1,bufr_nchan + scaleloop: do j=jstart,4 + if(allchan(1,i) >= cscale(1,j) .and. allchan(1,i) <= cscale(2,j))then + scalef(i) = sscale(j) + jstart=j + exit scaleloop + end if + end do scaleloop + end do + + if (iret /= bufr_nchan) then + write(6,*)'READ_IASI-NG: ### ERROR IN READING ', senname, ' BUFR DATA:', & + iret, ' CH DATA IS READ INSTEAD OF ',bufr_nchan + cycle read_loop + endif + +! Coordinate bufr channels with satinfo file channels +! If this is the first time or a change in the bufr channels is detected, sync with satinfo file + if (ANY(int(allchan(1,:)) /= bufr_chan_test(:))) then + sfc_channel_index = 0 + bufr_index(:) = 0 + bufr_chans: do l=1,bufr_nchan + bufr_chan_test(l) = int(allchan(1,l)) ! Copy this bufr channel selection into array for comparison to next profile + satinfo_chans: do i=1,satinfo_nchan ! Loop through sensor (iasi-ng) channels in the satinfo file + if ( channel_number(i) == int(allchan(1,l)) ) then ! Channel found in both bufr and satinfo file + bufr_index(i) = l + if ( channel_number(i) == sfc_channel) sfc_channel_index = l + exit satinfo_chans ! go to next bufr channel + endif + end do satinfo_chans + end do bufr_chans + endif + + if (sfc_channel_index == 0) then + write(6,*)'READ_IASI-NG: ***ERROR*** SURFACE CHANNEL USED FOR QC WAS NOT FOUND' + cycle read_loop + endif + +!$omp parallel do schedule(dynamic,1) private(i,sc_chan,bufr_chan,radiance) + channel_loop: do i=1,satinfo_nchan + bufr_chan = bufr_index(i) + if (bufr_chan > 0 ) then +! check that channel number is within reason + if (( allchan(2,bufr_chan) > zero .and. allchan(2,bufr_chan) < 99999._r_kind)) then ! radiance bounds + radiance = allchan(2,bufr_chan)*scalef(bufr_chan) + sc_chan = sc_index(i) + call crtm_planck_temperature(sensorindex_iasing,sc_chan,radiance,temperature(bufr_chan)) + else + temperature(bufr_chan) = tbmin + endif + end if + end do channel_loop + +! Check for reasonable temperature values + iskip = 0 + skip_loop: do i=1,satinfo_nchan + if ( bufr_index(i) == 0 ) cycle skip_loop + bufr_chan = bufr_index(i) + if(temperature(bufr_chan) <= tbmin .or. temperature(bufr_chan) > tbmax ) then + temperature(bufr_chan) = min(tbmax,max(tbmin,temperature(bufr_chan))) + if(iuse_rad(ioff+i) >= 0)iskip = iskip + 1 + endif + end do skip_loop + + if(iskip > 0)then + if(print_verbose)write(6,*) ' READ_IASI-NG: iskip > 0 ',iskip + cycle read_loop + end if + +! crit1=crit1 + ten*real(iskip,r_kind) + +! If the surface channel exists (~960.0 cm-1) and the imager cloud information is missing, use an +! estimate of the surface temperature to determine if the profile may be clear. + if (.not. cloud_info) then + pred = tsavg*0.98_r_kind - temperature(sfc_channel_index) + pred = max(pred,zero) + crit1=crit1 + pred + endif + +! Map obs to grids + if (pred == zero) then + call finalcheck(dist1,crit1,itx,iuse) + else + call finalcheck(one,crit1,itx,iuse) + endif + if(.not. iuse)cycle read_loop + +! Read the imager cluster information for the Cloud and Aerosol Detection Software. +! Only channels 4 and 5 are used. + + if ( iasing_cads ) then + call ufbseq(lnbufr,imager_info,123,7,iret,'IASICSSQ') + 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 + 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(:)) + imager_conversion(1) = one / (sc(sensorindex_imager)%wavenumber(18) **2) + imager_conversion(2) = one / (sc(sensorindex_imager)%wavenumber(19) **2) + +! Order clusters from largest (1) to smallest (7) + imager_cluster_sort: do i=1,7 + j = maxloc(imager_cluster_size,dim=1,mask=imager_cluster_flag) + imager_cluster_index(i) = j + imager_cluster_flag(j) = .FALSE. + end do imager_cluster_sort + +! Convert from radiance to brightness temperature for mean and standard deviation used by CADS. +! Imager cluster info added to data_all array + + 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) + + iexponent = -(nint(imager_info(104,i)) -11 ) ! channel 4 radiance for each cluster. + imager_info(105,i) = imager_info(105,i) * imager_conversion(1) * (ten ** iexponent) + + iexponent = -(nint(imager_info(106,i)) -11 ) ! channel 4 radiance std dev for each cluster. + imager_info(107,i) = imager_info(107,i) * imager_conversion(1) * (ten ** iexponent) + + call crtm_planck_temperature(sensorindex_imager,18,imager_info(105,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(110,i)) -11 ) ! channel 5 radiance for each cluster + imager_info(111,i) = imager_info(111,i) * imager_conversion(2) * (ten ** iexponent) + + iexponent = -(nint(imager_info(112,i)) -11 ) ! channel 5 radiance std dev for each cluser. + imager_info(113,i) = imager_info(113,i) * imager_conversion(2) * (ten ** iexponent) + + call crtm_planck_temperature(sensorindex_imager,19,imager_info(110,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 + +! Compute cluster averages for each channel + + imager_mean(1) = sum(imager_cluster_size(:) * imager_info(105,:)) ! Channel 4 radiance cluster average + imager_std_dev(1) = sum(imager_cluster_size(:) * (imager_info(105,:)**2 + imager_info(107,:)**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(111,:)) ! Channel 5 radiance cluster average + imager_std_dev(2) = sum(imager_cluster_size(:) * (imager_info(111,:)**2 + imager_info(113,:)**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) + + else ! Imager cluster information is missing. Set everything to zero + data_all(maxinfo+1 : maxinfo+25,itx) = zero + endif + endif ! iasing_cads = .true. +! +! interpolate NSST variables to Obs. location and get dtw, dtc, tz_tr +! + if ( nst_gsi > 0 ) then + tref = ts(0) + dtw = zero + dtc = zero + tz_tr = one + if ( sfcpct(0) > zero ) then + call gsi_nstcoupler_deter(dlat_earth,dlon_earth,t4dv,zob,tref,dtw,dtc,tz_tr) + endif + endif + + rsat=allspot(1) + data_all(1,itx) = rsat ! satellite ID + data_all(2,itx) = t4dv ! time diff (obs-anal) (hrs) + data_all(3,itx) = dlon ! grid relative longitude + data_all(4,itx) = dlat ! grid relative latitude + data_all(5,itx) = sat_zenang*deg2rad ! satellite zenith angle (rad) + data_all(6,itx) = allspot(11) ! satellite azimuth angle (deg) + data_all(7,itx) = lza ! look angle (rad) + data_all(8,itx) = ifovn ! fov number + data_all(9,itx) = allspot(12) ! solar zenith angle (deg) + data_all(10,itx)= allspot(13) ! solar azimuth angle (deg) + data_all(11,itx) = sfcpct(0) ! sea percentage of + data_all(12,itx) = sfcpct(1) ! land percentage + data_all(13,itx) = sfcpct(2) ! sea ice percentage + data_all(14,itx) = sfcpct(3) ! snow percentage + data_all(15,itx)= ts(0) ! ocean skin temperature + data_all(16,itx)= ts(1) ! land skin temperature + data_all(17,itx)= ts(2) ! ice skin temperature + data_all(18,itx)= ts(3) ! snow skin temperature + data_all(19,itx)= tsavg ! average skin temperature + data_all(20,itx)= vty ! vegetation type + data_all(21,itx)= vfr ! vegetation fraction + data_all(22,itx)= sty ! soil type + data_all(23,itx)= stp ! soil temperature + data_all(24,itx)= sm ! soil moisture + data_all(25,itx)= sn ! snow depth + data_all(26,itx)= zz ! surface height + data_all(27,itx)= idomsfc(1) + 0.001_r_kind ! dominate surface type + data_all(28,itx)= sfcr ! surface roughness + data_all(29,itx)= ff10 ! ten meter wind factor + data_all(30,itx)= dlon_earth_deg ! earth relative longitude (degrees) + data_all(31,itx)= dlat_earth_deg ! earth relative latitude (degrees) + + if(dval_use)then + data_all(maxinfo+cads_info+1,itx)= val_iasing + data_all(maxinfo+cads_info+2,itx)= itt + end if + + if ( nst_gsi > 0 ) then + data_all(maxinfo+cads_info+dval_info+1,itx) = tref ! foundation temperature + data_all(maxinfo+cads_info+dval_info+2,itx) = dtw ! dt_warm at zob + data_all(maxinfo+cads_info+dval_info+3,itx) = dtc ! dt_cool at zob + data_all(maxinfo+cads_info+dval_info+4,itx) = tz_tr ! d(Tz)/d(Tr) + endif + +! Put satinfo defined channel temperatures into data array + do l=1,satinfo_nchan + ! Prevent out of bounds reference from temperature + if ( bufr_index(l) == 0 ) cycle + i = bufr_index(l) + if(i /= 0)then + data_all(l+nreal,itx) = temperature(i) ! brightness temerature + else + data_all(l+nreal,itx) = tbmin + end if + end do + nrec(itx)=irec + + enddo read_loop + + enddo read_subset + + call closbf(lnbufr) + close(lnbufr) + + end do ears_db_loop + + deallocate(temperature, allchan, bufr_chan_test,scalef) + deallocate(channel_number,sc_index) + deallocate(bufr_index) +! deallocate crtm info + error_status = crtm_spccoeff_destroy() + if (error_status /= success) & + write(6,*)'OBSERVER: ***ERROR*** crtm_destroy error_status=',error_status + +! If multiple tasks read input bufr file, allow each tasks to write out +! information it retained and then let single task merge files together + + call combine_radobs(mype_sub,mype_root,npe_sub,mpi_comm_sub,& + nele,itxmax,nread,ndata,data_all,score_crit,nrec) + +! Allow single task to check for bad obs, update superobs sum, +! and write out data to scratch file for further processing. + if (mype_sub==mype_root.and.ndata>0) then + +! Identify "bad" observation (unreasonable brightness temperatures). +! Update superobs sum according to observation location + + do n=1,ndata + do i=1,satinfo_nchan + if(data_all(i+nreal,n) > tbmin .and. & + data_all(i+nreal,n) < tbmax)nodata=nodata+1 + end do + end do + + if(dval_use .and. assim)then + do n=1,ndata + itt=nint(data_all(33,n)) + super_val(itt)=super_val(itt)+val_iasing + end do + end if + +! Write final set of "best" observations to output file + call count_obs(ndata,nele,ilat,ilon,data_all,nobs) + write(lunout) obstype,sis,nreal,satinfo_nchan,ilat,ilon + write(lunout) ((data_all(k,n),k=1,nele),n=1,ndata) + + endif + + + deallocate(data_all,nrec) ! Deallocate data arrays + call destroygrids ! Deallocate satthin arrays + +! Deallocate arrays and nullify pointers. + if(isfcalc == 1) call fov_cleanup + + if(diagnostic_reg .and. ntest > 0 .and. mype_sub==mype_root) & + write(6,*)'READ_IASI-NG: mype,ntest,disterrmax=',& + mype,ntest,disterrmax + + return +end subroutine read_iasing diff --git a/src/gsi/read_obs.F90 b/src/gsi/read_obs.F90 index aa0f11b4e3..17e8db9be4 100644 --- a/src/gsi/read_obs.F90 +++ b/src/gsi/read_obs.F90 @@ -920,7 +920,7 @@ subroutine read_obs(ndata,mype) ditype(i) = 'wcp' else if( hirs .or. sndr .or. seviri .or. abi .or. & obstype == 'airs' .or. obstype == 'amsua' .or. & - obstype == 'msu' .or. obstype == 'iasi' .or. & + obstype == 'msu' .or. obstype == 'iasi' .or. obstype == 'iasi-ng' .or. & obstype == 'amsub' .or. obstype == 'mhs' .or. & obstype == 'hsb' .or. obstype == 'goes_img' .or. & obstype == 'ahi' .or. avhrr .or. & @@ -1013,6 +1013,8 @@ subroutine read_obs(ndata,mype) parallel_read(i)= .true. else if(obstype == 'iasi')then parallel_read(i)= .true. + else if(obstype == 'iasi-ng')then + parallel_read(i)= .true. else if(obstype == 'amsub')then parallel_read(i)= .true. else if(obstype == 'mhs' )then @@ -1062,23 +1064,24 @@ subroutine read_obs(ndata,mype) (obstype == 'amsua' .or. obstype == 'amsub' .or. & obstype == 'mhs' .or. obstype == 'hirs3' .or. & obstype == 'cris' .or. obstype == 'cris-fsr' .or. & - obstype == 'iasi' .or. obstype == 'atms') .and. & + obstype == 'iasi' .or. obstype == 'iasi-ng' .or. & + obstype == 'atms') .and. & (dplat(i) == 'n17' .or. dplat(i) == 'n18' .or. & dplat(i) == 'n19' .or. dplat(i) == 'npp' .or. & dplat(i) == 'n20' .or. dplat(i) == 'n21' .or. & dplat(i) == 'metop-a' .or. dplat(i) == 'metop-b' .or. & - dplat(i) == 'metop-c') + dplat(i) == 'metop-c' .or. dplat(i) == 'metop-sg-a1') ! direct broadcast from NESDIS/UW db_possible(i) = ditype(i) == 'rad' .and. & (obstype == 'amsua' .or. obstype == 'amsub' .or. & obstype == 'mhs' .or. obstype == 'atms' .or. & obstype == 'cris' .or. obstype == 'cris-fsr' .or. & - obstype == 'iasi') .and. & + obstype == 'iasi' .or. obstype == 'iasi-ng') .and. & (dplat(i) == 'n17' .or. dplat(i) == 'n18' .or. & dplat(i) == 'n19' .or. dplat(i) == 'npp' .or. & dplat(i) == 'n20' .or. dplat(i) == 'n21' .or. & dplat(i) == 'metop-a' .or. dplat(i) == 'metop-b' .or. & - dplat(i) == 'metop-c') + dplat(i) == 'metop-c' .or. dplat(i) == 'metop-sg-a1') ! Inquire data set to deterimine if input data available and size of dataset ii=ii+1 @@ -1743,6 +1746,14 @@ subroutine read_obs(ndata,mype) read_rec(i),read_ears_rec(i),read_db_rec(i),dval_use) string='READ_IASI' +! Process iasi-ng data + else if(obstype == 'iasi-ng')then + call read_iasing(mype,val_dat,ithin,isfcalc,rmesh,platid,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) + string='READ_IASI-NG' + ! Process cris data else if(obstype == 'cris' .or. obstype =='cris-fsr' )then call read_cris(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,& diff --git a/src/gsi/setuprad.f90 b/src/gsi/setuprad.f90 index 136568d1f3..7d14f762f8 100644 --- a/src/gsi/setuprad.f90 +++ b/src/gsi/setuprad.f90 @@ -297,7 +297,7 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& use radinfo, only: iland_det, isnow_det, iwater_det, imix_det, iice_det, & iomg_det, itopo_det, isst_det,iwndspeed_det, optconv use qcmod, only: setup_tzr_qc,ifail_scanedge_qc,ifail_outside_range - use qcmod, only: iasi_cads, cris_cads + use qcmod, only: iasi_cads, iasing_cads, cris_cads use state_vectors, only: svars3d, levels, svars2d, ns3d use oneobmod, only: lsingleradob,obchan,oblat,oblon,oneob_type use correlated_obsmod, only: corr_adjust_jacobian, idnames @@ -367,7 +367,7 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& logical cao_flag logical hirs2,msu,goessndr,hirs3,hirs4,hirs,amsua,amsub,airs,hsb,goes_img,ahi,mhs,abi type(sparr2) :: dhx_dx - logical avhrr,avhrr_navy,viirs,lextra,ssu,iasi,cris,seviri,atms + logical avhrr,avhrr_navy,viirs,lextra,ssu,iasi,iasing,cris,seviri,atms logical ssmi,ssmis,amsre,amsre_low,amsre_mid,amsre_hig,amsr2,gmi,saphir logical ssmis_las,ssmis_uas,ssmis_env,ssmis_img logical sea,mixed,land,ice,snow,toss,l_may_be_passive,eff_area @@ -522,6 +522,7 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& ssmis_img = obstype == 'ssmis_img' ssmis_env = obstype == 'ssmis_env' iasi = obstype == 'iasi' + iasing = obstype == 'iasi-ng' cris = obstype == 'cris' .or. obstype == 'cris-fsr' seviri = obstype == 'seviri' atms = obstype == 'atms' @@ -610,7 +611,7 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& imager_cluster_bt=zero imager_chan_stdev=zero imager_model_bt=zero - if ((iasi_cads .and. iasi) .or. (cris_cads .and. cris)) then + if ((iasi_cads .and. iasi) .or. (iasing_cads .and. iasing) .or. (cris_cads .and. cris)) then call cads_imager_calc(obstype,isis,nobs,nreal,nchanl,nsig,data_s,init_pass,mype, & imager_cluster_fraction,imager_cluster_bt,imager_chan_stdev, imager_model_bt) @@ -1355,7 +1356,7 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& ! ---------- IR ------------------- ! QC HIRS/2, GOES, HIRS/3 and AIRS sounder data ! - ObsQCs: if (hirs .or. goessndr .or. airs .or. iasi .or. cris) then + ObsQCs: if (hirs .or. goessndr .or. airs .or. iasi .or. iasing .or. cris) then frac_sea=data_s(ifrac_sea,n) @@ -1376,7 +1377,7 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& end do call qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse(n),goessndr,airs,cris,iasi, & - hirs,zsges,cenlat,frac_sea,pangs,trop5,zasat,tzbgr,tsavg5,tbc,tb_obs,tbcnob,tnoise, & + iasing,hirs,zsges,cenlat,frac_sea,pangs,trop5,zasat,tzbgr,tsavg5,tbc,tb_obs,tbcnob,tnoise, & wavenumber,ptau5,prsltmp,tvp,temp,wmix,chan_level,emissivity_k,ts,tsim, & id_qc,aivals,errf,varinv,varinv_use,cld,cldp,kmax,zero_irjaco3_pole(n), & imager_cluster_fraction(:,n), imager_cluster_bt(:,:,n), imager_chan_stdev(:,n),imager_model_bt(:,n)) @@ -1791,7 +1792,7 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& varinv0(ii)=varinv(ii) raterr2(ii)=error0(ii)**2*varinv0(ii) if (l_may_be_passive .and. .not. retrieval) then - if(optconv > zero .and. (iasi .or. cris) .and. iinstr /= -1)then + if(optconv > zero .and. (iasi .or. iasing .or. cris) .and. iinstr /= -1)then asum=zero do k=1,nsig asum=asum+abs(jacobian(iqs+k,ii))*qs(k) From 0e1853b9fd8b3dd57820e7760f211f9b138510fd Mon Sep 17 00:00:00 2001 From: wx20jjung Date: Wed, 10 Jul 2024 23:33:43 +0000 Subject: [PATCH 03/16] First save for iasi-ng. added read_iasing and various tests for iasi-ng --- src/gsi/crtm_interface.f90 | 3 +- src/gsi/gsi_obOperTypeManager.F90 | 1 + src/gsi/gsimod.F90 | 2 +- src/gsi/mrmsmod.f90 | 5 +- src/gsi/obsmod.F90 | 3 +- src/gsi/qcmod.f90 | 9 +-- src/gsi/radiance_mod.f90 | 3 +- src/gsi/radinfo.f90 | 11 ++-- src/gsi/read_cris.f90 | 6 +- src/gsi/read_diag.f90 | 8 ++- src/gsi/read_iasi.f90 | 2 +- src/gsi/read_iasing.f90 | 96 ++++++++++++++++++++----------- src/gsi/read_obs.F90 | 63 ++++++++++---------- src/gsi/setuprad.f90 | 2 +- src/gsi/statsrad.f90 | 6 +- 15 files changed, 129 insertions(+), 91 deletions(-) diff --git a/src/gsi/crtm_interface.f90 b/src/gsi/crtm_interface.f90 index 8a7fb84507..192070463e 100644 --- a/src/gsi/crtm_interface.f90 +++ b/src/gsi/crtm_interface.f90 @@ -662,7 +662,8 @@ subroutine init_crtm(init_pass,mype_diaghdr,mype,nchanl,nreal,isis,obstype,radmo error_status = crtm_channelinfo_subset(channelinfo(1), & channel_subset = nuchan(subset_start:subset_end)) -else if (channelinfo(1)%sensor_id(1:7) == 'iasi-ng' .AND. isis(1:7) == 'iasi-ng') then +!JAJ the correct one else if (channelinfo(1)%sensor_id(1:7) == 'iasi-ng' .AND. isis(1:7) == 'iasi-ng') then +else if (channelinfo(1)%sensor_id(1:3) == '999' .AND. isis(1:7) == 'iasi-ng') then sensorindex = 1 subset_start = 0 subset_end = 0 diff --git a/src/gsi/gsi_obOperTypeManager.F90 b/src/gsi/gsi_obOperTypeManager.F90 index 6db7921905..3f6dce3a0a 100644 --- a/src/gsi/gsi_obOperTypeManager.F90 +++ b/src/gsi/gsi_obOperTypeManager.F90 @@ -321,6 +321,7 @@ function dtype2index_(dtype) result(index_) case("hsb" ); index_= iobOper_rad ! case("iasi" ); index_= iobOper_rad + case("iasi-ng"); index_= iobOper_rad case("cris" ); index_= iobOper_rad case("cris-fsr" ); index_= iobOper_rad ! diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index d11616b55d..99bedda7e9 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -2321,7 +2321,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,a12,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..e72158baf2 100644 --- a/src/gsi/mrmsmod.f90 +++ b/src/gsi/mrmsmod.f90 @@ -27,7 +27,7 @@ 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, ditype, dplat character(20),intent(inout),dimension(nrows):: obsfile_all character(*),intent(inout),dimension(nrows):: dfile character(20),intent(inout),dimension(nrows):: dsis @@ -42,7 +42,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 95e5046d9f..e73064be90 100644 --- a/src/gsi/obsmod.F90 +++ b/src/gsi/obsmod.F90 @@ -625,7 +625,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 6425d5268e..f1e81323c8 100644 --- a/src/gsi/qcmod.f90 +++ b/src/gsi/qcmod.f90 @@ -2335,9 +2335,9 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr,airs boundary_layer_pres = nint(0.8_r_kind*prsltmp(1)) ! boundary layer set to be 80% of surface pressure tropopause_height = nint(trop5) imager_chans = (/18,19/) ! imager channel numbers (from satinfo) - isurface_chan = 2447 ! surface channel - ichan_10_micron = 2447 ! ~10.7 micron channel for low level cloud test - ichan_12_micron = 1560 ! ~12.0 micron channel for low level cloud test + isurface_chan = 2539 ! surface channel + ichan_10_micron = 2343 ! ~10.7 micron channel for low level cloud test + ichan_12_micron = 1509 ! ~12.0 micron channel for low level cloud test call cloud_aerosol_detection( I_Sensor_ID, nchanl, chan_array, & tropopause_height, boundary_layer_pres, tb_bc, tsim, chan_level, imager_chans, cluster_fraction, & @@ -2434,6 +2434,7 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr,airs ! default compute cloud stats, emc_legacy_cloud_detect else if ( lcloud > 0 ) then +!JAJ if ( lcloud > 0 .and. .not. iasing ) then do i=1,nchanl ! reject channels with iuse_rad(j)=-1 when they are peaking below the cloud @@ -4258,7 +4259,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/radiance_mod.f90 b/src/gsi/radiance_mod.f90 index aae7794957..ee57a294a0 100644 --- a/src/gsi/radiance_mod.f90 +++ b/src/gsi/radiance_mod.f90 @@ -412,7 +412,8 @@ subroutine radiance_obstype_init rtype(i) == 'avhrr' .or. rtype(i) == 'amsre' .or. rtype(i) == 'ssmis' .or. & rtype(i) == 'ssmi' .or. rtype(i) == 'atms' .or. rtype(i) == 'cris' .or. & rtype(i) == 'amsr2' .or. rtype(i) == 'gmi' .or. rtype(i) == 'saphir' .or. & - rtype(i) == 'cris-fsr' .or. rtype(i) == 'abi' .or. rtype(i) == 'viirs' ) then + rtype(i) == 'cris-fsr' .or. rtype(i) == 'abi' .or. rtype(i) == 'viirs' .or. & + rtype(i) == 'iasi-ng' )then drtype(i)='rads' end if end do diff --git a/src/gsi/radinfo.f90 b/src/gsi/radinfo.f90 index 4ad17626e6..e47ee7def4 100644 --- a/src/gsi/radinfo.f90 +++ b/src/gsi/radinfo.f90 @@ -768,7 +768,8 @@ subroutine radinfo_read end select if ( .not. diag_rad .and. iuse_rad(j) < 0 .and. iextra_det(j) < 0 .and. & - ( nusis(j)(1:4) == 'cris' .or. nusis(j)(1:4) == 'iasi' .or. nusis(j)(1:4) == 'airs')) cycle + ( nusis(j)(1:4) == 'cris' .or. nusis(j)(1:7) == 'iasi-ng' .or. nusis(j)(1:4) == 'iasi' & + .or. nusis(j)(1:4) == 'airs')) cycle if(iuse_rad(j) == 4 .or. iuse_rad(j) == 2) air_rad(j)=zero if(iuse_rad(j) == 4 .or. iuse_rad(j) == 3) ang_rad(j)=zero @@ -805,10 +806,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) @@ -869,7 +870,7 @@ subroutine radinfo_read ! The second part of the if statement keeps from printing them. if ( .not. cfound ) then if ((diag_rad .and. mype ==0) .or. & - (.not. diag_rad .and. isis(1:4)/='airs' .and. isis(1:4) /= 'cris' .and. isis(1:4) /= 'iasi')) & + (.not. diag_rad .and. isis(1:4)/='airs' .and. isis(1:4) /= 'cris' .and. isis(1:7) /= 'iasi-ng' .and. isis(1:4) /= 'iasi')) & write(6,*) '***WARNING instrument/channel ',isis,ichan,'found in satbias_pc file but not found in satinfo' endif @@ -1112,7 +1113,7 @@ subroutine radinfo_read ! The second part of the if statement keeps from printing them. if ( .not. cfound ) then if ((diag_rad .and. mype ==0) .or. & - (.not. diag_rad .and. isis(1:4)/='airs' .and. isis(1:4) /= 'cris' .and. isis(1:4) /= 'iasi')) & + (.not. diag_rad .and. isis(1:4)/='airs' .and. isis(1:4) /= 'cris' .and. isis(1:7) /= 'iasi-ng' .and. isis(1:4) /= 'iasi')) & write(6,*) '***WARNING instrument/channel ',isis,ichan,'found in satbias_in file but not found in satinfo' endif diff --git a/src/gsi/read_cris.f90 b/src/gsi/read_cris.f90 index 7e9195fe41..3341dfbed9 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 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..e8ef9a2a58 100644 --- a/src/gsi/read_iasi.f90 +++ b/src/gsi/read_iasi.f90 @@ -143,7 +143,7 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& 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 ) :: jsatid character(len=*), intent(in ) :: obstype character(len=20),intent(in ) :: sis real(r_kind) ,intent(in ) :: twind diff --git a/src/gsi/read_iasing.f90 b/src/gsi/read_iasing.f90 index d0ebb9d9c1..af52d36456 100644 --- a/src/gsi/read_iasing.f90 +++ b/src/gsi/read_iasing.f90 @@ -93,9 +93,9 @@ subroutine read_iasing(mype,val_iasing,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, obstype, jsatid +! 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_iasing @@ -160,7 +160,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& logical :: outside,iuse,assim,valid logical :: iasing,quiet,cloud_info - integer(i_kind) :: ifov, instr, iscn, ioff, sensorindex_iasing + integer(i_kind) :: ifov, ifor, istep, ipos, instr, iscn, ioff, sensorindex_iasing integer(i_kind) :: i, j, l, iskip, ifovn, bad_line, ksatid, kidsat, llll integer(i_kind) :: nreal, isflg integer(i_kind) :: itx, k, nele, itt, n @@ -188,7 +188,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& ! Set standard parameters character(8),parameter:: fov_flag="crosstrk" - integer(i_kind),parameter:: sfc_channel=1271 + integer(i_kind),parameter:: sfc_channel=2539 integer(i_kind),parameter:: ichan=-999 ! fov-based surface code is not channel specific for iasing real(r_kind),parameter:: expansion=one ! exansion factor for fov-based surface code. ! use one for ir sensors. @@ -226,7 +226,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& call gsi_nstcoupler_skindepth(obstype, zob) ! get penetration depth (zob) for the obstype endif - if(jsatid == 'metop-sg-a1')kidsat=24 + if (jsatid == 'metop-sg-a1') kidsat=24 ! write(6,*)'READ_IASI-NG: mype, mype_root,mype_sub, npe_sub,mpi_comm_sub', & ! mype, mype_root,mype_sub,mpi_comm_sub @@ -244,8 +244,8 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& if (trim(nusis(i))==trim(sis)) then ioff = min(ioff,i) ! iasi-ng offset if (subset_start == 0) then - step = radstep(i) - start = radstart(i) +!JAJ step = radstep(i) +!JAJ start = radstart(i) if (radedge1(i)/=-1 .and. radedge2(i)/=-1) then radedge_min=radedge1(i) radedge_max=radedge2(i) @@ -307,7 +307,8 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& ! find IASI-NG sensorindex sensorindex_iasing = 0 - if ( sc(1)%sensor_id(1:7) == 'iasi-ng' ) then +!JAJ if ( sc(1)%sensor_id(1:7) == 'iasi-ng' ) then + if ( sc(1)%sensor_id == '999' ) then sensorindex_iasing = 1 else write(6,*)'READ_IASI-NG: ***ERROR*** sensorindex_iasi-ng not set NO IASI-NG DATA USED' @@ -317,6 +318,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& ! find imager sensorindex sensorindex_imager = 0 + iasing_cads = .true. if ( iasing_cads .and. imager_coeff ) then if ( sc(2)%sensor_id(1:7) == 'metimag' ) then sensorindex_imager = 2 @@ -459,20 +461,20 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& ifov = nint(linele(1)) ! field of view -! IASI-NG fov ranges from 1 to 120. Current angle dependent bias +! IASI-NG fov ranges from 1 to 112. Current angle dependent bias ! correction has a maximum of 90 scan positions. Geometry -! of IASI-NG scan allows us to remap 1-120 to 1-60. Variable +! of IASI-NG scan allows us to remap 1-112 to 1-56. Variable ! ifovn below contains the remapped IASI-NG fov. This value is ! passed on to and used in setuprad - ifovn = (ifov-1)/2 + 1 +!JAJ ifovn = (ifov-1)/2 + 1 ! Remove data on edges - if (.not. use_edges .and. & - (ifovn < radedge_min .OR. ifovn > radedge_max )) cycle read_loop +!JAJ if (.not. use_edges .and. & +!JAJ (ifovn < radedge_min .OR. ifovn > radedge_max )) cycle read_loop ! Check field of view (FOVN) and satellite zenith angle (SAZA) iscn = nint(linele(2)) ! scan line - if( ifov <= 0 .or. ifov > 120) then + if( ifov <= 0 .or. ifov > 16) then write(6,*)'READ_IASI-NG: ### ERROR IN READING ', senname, ' BUFR DATA:', & ' STRANGE OBS INFO(FOVN,SLNM):', ifov, iscn cycle read_loop @@ -481,7 +483,6 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& call ufbint(lnbufr,allspot,13,1,iret,allspotlist) if(iret /= 1) cycle read_loop - ! Check observing position dlat_earth = allspot(8) ! latitude dlon_earth = allspot(9) ! longitude @@ -583,21 +584,39 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& ' STRANGE OBS INFO(FOVN,SLNM,SAZA,BEARAZ):', ifov, iscn, allspot(10),allspot(11) cycle read_loop endif - if ( ifov <= 60 ) sat_zenang = -sat_zenang - -! Compare IASI-NG satellite scan angle and zenith angle - piece = -step_adjust - if ( mod(ifovn,2) == 1) piece = step_adjust - lza = ((start + real((ifov-1)/4,r_kind)*step) + piece)*deg2rad - sat_height_ratio = (earth_radius + linele(4))/earth_radius - lzaest = asin(sat_height_ratio*sin(lza))*rad2deg - if (abs(sat_zenang - lzaest) > one) then - write(6,*)' READ_IASI-NG WARNING uncertainty in lza ', & - lza*rad2deg,sat_zenang,sis,ifov,start,step,allspot(11),allspot(12),allspot(13) - bad_line = iscn - cycle read_loop - endif - +! Remove this code when the field-of-view spans 1-224 JAJ. + ifov = ( 59 - nint(sat_zenang)) / 2 + ifov = min( 28, max(1,ifov)) + +! There are 16 fields-of-view within each field-of-reagard, a 4X4 matrix. +! | 1 2 3 4 | +! | 5 6 7 8 | +! | 9 10 11 12 | +! | 13 14 15 16 | + +! There are 14 fields-of-regard in each scan line. So, there are 56 unique positions in the scan line. +! To determine the scan position: + ifor = ((ifov-1) / 16) + 1 ! Determine field-of-regard + istep = ifov - (ifor * 16) ! Determine field-of-view within field-of-regard + ipos = (ifor * 4) + mod(istep,4) + 1 ! Determine position of field-of-view within scan line + +! Remove data on edges + if (.not. use_edges .and. & + (ipos < radedge_min .or. ipos > radedge_max)) cycle read_loop + + if ( ipos <= 28 ) sat_zenang = - sat_zenang + +!JAJ to be removed +! Estimate scan angle and compare to zenith angle +!JAJ lza = ((start + real(ipos)* step) * deg2rad +!JAJ sat_height_ratio = (earth_radius + linele(4)) / earth_radius +!JAJ lzaest = asin(sat_height_ratio * sin(lza)) * rad2deg +!JAJ if (abs(sat_zenang - lzaest) > one) then +!JAJ write(6,*)' READ_IASI-NG WARNING uncertainty in lza ', & +!JAJ lza*rad2deg,sat_zenang,sis,ifov,start,step,allspot(11),allspot(12),allspot(13) +!JAJ bad_line = iscn +!JAJ cycle read_loop +!JAJ endif ! "Score" observation. We use this information to identify "best" obs ! Locate the observation on the analysis grid. Get sst and land/sea/ice @@ -708,7 +727,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& cycle read_loop endif -!$omp parallel do schedule(dynamic,1) private(i,sc_chan,bufr_chan,radiance) +!JAJ !$omp parallel do schedule(dynamic,1) private(i,sc_chan,bufr_chan,radiance) channel_loop: do i=1,satinfo_nchan bufr_chan = bufr_index(i) if (bufr_chan > 0 ) then @@ -717,6 +736,12 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& radiance = allchan(2,bufr_chan)*scalef(bufr_chan) sc_chan = sc_index(i) call crtm_planck_temperature(sensorindex_iasing,sc_chan,radiance,temperature(bufr_chan)) +!JAJ if (( temperature(bufr_chan) > 220.0_r_kind .and. temperature(bufr_chan) < 300.0_r_kind) .and. & +!JAJ (bufr_chan > 500 .and. bufr_chan < 5000)) then +!JAJ write(*,'(a10,1x,i6,1x,e15.3,1x,f7.2,1x,f10.5,1x,f10.5)') 'JAJ value ',bufr_chan,allchan(2,bufr_chan),temperature(bufr_chan), & +!JAJ dlat_earth_deg, dlon_earth_deg +!JAJ endif + else temperature(bufr_chan) = tbmin endif @@ -730,6 +755,8 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& bufr_chan = bufr_index(i) if(temperature(bufr_chan) <= tbmin .or. temperature(bufr_chan) > tbmax ) then temperature(bufr_chan) = min(tbmax,max(tbmin,temperature(bufr_chan))) + !write(*,'(a10,1x,i6,1x,e7.3,1x,f9.5,1x,f9.5)') 'JAJ bad value ',bufr_chan,allchan(2,bufr_chan),scalef(bufr_chan) +! write(*,'(a14,1x,i6,1x,e12.3,1x,f10.5,1x,f10.5)') 'JAJ bad value ',bufr_chan,allchan(2,bufr_chan),dlat_earth_deg, dlon_earth_deg if(iuse_rad(ioff+i) >= 0)iskip = iskip + 1 endif end do skip_loop @@ -851,8 +878,9 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& data_all(4,itx) = dlat ! grid relative latitude data_all(5,itx) = sat_zenang*deg2rad ! satellite zenith angle (rad) data_all(6,itx) = allspot(11) ! satellite azimuth angle (deg) - data_all(7,itx) = lza ! look angle (rad) - data_all(8,itx) = ifovn ! fov number +!JAJ data_all(7,itx) = lza ! look angle (rad) + data_all(7,itx) = sat_zenang*deg2rad ! look angle (rad) + data_all(8,itx) = ifov ! fov number data_all(9,itx) = allspot(12) ! solar zenith angle (deg) data_all(10,itx)= allspot(13) ! solar azimuth angle (deg) data_all(11,itx) = sfcpct(0) ! sea percentage of diff --git a/src/gsi/read_obs.F90 b/src/gsi/read_obs.F90 index 17e8db9be4..59d736e067 100644 --- a/src/gsi/read_obs.F90 +++ b/src/gsi/read_obs.F90 @@ -229,6 +229,8 @@ subroutine read_obs_check (lexist,filename,jsatid,dtype,minuse,nread) kidsat=3 else if(jsatid == 'metop-c')then kidsat=5 + else if(jsatid == 'metop-sg-a1')then + kidsat=24 else if(jsatid == 'm08')then kidsat = 55 else if(jsatid == 'm09')then @@ -778,7 +780,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 @@ -1405,7 +1407,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 @@ -1538,12 +1539,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' @@ -1701,12 +1702,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) @@ -1714,7 +1715,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) @@ -1722,7 +1723,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) @@ -1730,9 +1731,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) @@ -1740,7 +1741,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) @@ -1748,7 +1749,7 @@ subroutine read_obs(ndata,mype) ! Process iasi-ng data else if(obstype == 'iasi-ng')then - call read_iasing(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,& + call read_iasing(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) @@ -1756,7 +1757,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) @@ -1767,7 +1768,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) @@ -1775,7 +1776,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) @@ -1784,7 +1785,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) @@ -1794,7 +1795,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) @@ -1802,7 +1803,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)) @@ -1810,7 +1811,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) @@ -1818,7 +1819,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) @@ -1826,14 +1827,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) @@ -1841,7 +1842,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) @@ -1850,7 +1851,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) @@ -1858,14 +1859,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) @@ -1878,7 +1879,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' @@ -1888,7 +1889,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 @@ -1914,7 +1915,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' @@ -1940,7 +1941,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(i) =',trim(dplat(i))) call warn('read_obs',' string =',trim(string)) endif diff --git a/src/gsi/setuprad.f90 b/src/gsi/setuprad.f90 index 7d14f762f8..8ed14b6fee 100644 --- a/src/gsi/setuprad.f90 +++ b/src/gsi/setuprad.f90 @@ -522,7 +522,7 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& ssmis_img = obstype == 'ssmis_img' ssmis_env = obstype == 'ssmis_env' iasi = obstype == 'iasi' - iasing = obstype == 'iasi-ng' + iasing = obstype == 'iasi-ng' cris = obstype == 'cris' .or. obstype == 'cris-fsr' seviri = obstype == 'seviri' atms = obstype == 'atms' diff --git a/src/gsi/statsrad.f90 b/src/gsi/statsrad.f90 index d42a53f7d6..aac671b7b8 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,2i6,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) From 8f347ffc1c3be26c86d356ab573c420602f21e9a Mon Sep 17 00:00:00 2001 From: wx20jjung Date: Wed, 31 Jul 2024 17:08:11 +0000 Subject: [PATCH 04/16] The cloud and aerosol detection software (CADS) uses imager information. In this case it is METImage. Specifically, the CRTM coefficent files are needed by CADS. The CRTM coefficient files used are determined by the satinfo file. The METImage entry in the satinfo file generates error messages in various parts of the GSI. These changes silence the error messages. At this time, there is no METImage assimilation code in the GSI. --- src/gsi/gsi_obOperTypeManager.F90 | 1 + src/gsi/read_obs.F90 | 4 +++- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/gsi/gsi_obOperTypeManager.F90 b/src/gsi/gsi_obOperTypeManager.F90 index 3f6dce3a0a..f1549364c1 100644 --- a/src/gsi/gsi_obOperTypeManager.F90 +++ b/src/gsi/gsi_obOperTypeManager.F90 @@ -352,6 +352,7 @@ function dtype2index_(dtype) result(index_) ! case("avhrr_navy"); index_= iobOper_rad case("avhrr" ); index_= iobOper_rad + case("metimage" ); index_= iobOper_rad case("viirs-m" ); index_= iobOper_rad case("tcp" ,"[tcpoper]" ); index_= iobOper_tcp diff --git a/src/gsi/read_obs.F90 b/src/gsi/read_obs.F90 index 59d736e067..dd6f56fb34 100644 --- a/src/gsi/read_obs.F90 +++ b/src/gsi/read_obs.F90 @@ -929,7 +929,7 @@ subroutine read_obs(ndata,mype) amsre .or. ssmis .or. obstype == 'ssmi' .or. & obstype == 'ssu' .or. obstype == 'atms' .or. & obstype == 'cris' .or. obstype == 'cris-fsr' .or. & - obstype == 'amsr2' .or. obstype == 'viirs-m' .or. & + obstype == 'amsr2' .or. obstype == 'viirs-m' .or. obstype == 'metimage' .or. & obstype == 'gmi' .or. obstype == 'saphir' ) then ditype(i) = 'rad' else if (is_extOzone(dfile(i),obstype,dplat(i))) then @@ -1038,6 +1038,8 @@ subroutine read_obs(ndata,mype) parallel_read(i)= .true. else if(obstype == 'viirs-m' )then parallel_read(i)= .true. + else if(obstype == 'metimage' )then + parallel_read(i)= .true. else if(amsre)then parallel_read(i)= .true. else if(obstype == 'goes_img' )then From edc4eff44e666ad70cb8e72957b64bb4d119dfef Mon Sep 17 00:00:00 2001 From: wx20jjung Date: Tue, 3 Sep 2024 14:43:55 +0000 Subject: [PATCH 05/16] A few more changes for IASI-NG --- src/gsi/calc_fov_crosstrk.f90 | 2 + src/gsi/read_iasing.f90 | 83 ++++++++++++++++++++--------------- 2 files changed, 49 insertions(+), 36 deletions(-) diff --git a/src/gsi/calc_fov_crosstrk.f90 b/src/gsi/calc_fov_crosstrk.f90 index 6cb817b56b..816e820890 100644 --- a/src/gsi/calc_fov_crosstrk.f90 +++ b/src/gsi/calc_fov_crosstrk.f90 @@ -1289,6 +1289,8 @@ subroutine get_sat_height(satid, height, valid) height=840._r_kind case('n20', 'n21', 'n22', 'n23') height=840._r_kind + case('metop-sg-a1', 'metop-sg-a2', 'metop-sg-a3') + height=840._r_kind case default write(6,*) 'GET_SAT_HEIGHT: ERROR, unrecognized satellite id: ', trim(satid) valid=.false. diff --git a/src/gsi/read_iasing.f90 b/src/gsi/read_iasing.f90 index af52d36456..ad3611de21 100644 --- a/src/gsi/read_iasing.f90 +++ b/src/gsi/read_iasing.f90 @@ -94,8 +94,6 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& integer(i_kind) ,intent(in ) :: npe_sub integer(i_kind) ,intent(in ) :: mpi_comm_sub character(len=*), intent(in ) :: infile, obstype, jsatid -! 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_iasing @@ -184,7 +182,6 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& real(r_kind),dimension(123,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 ! Set standard parameters character(8),parameter:: fov_flag="crosstrk" @@ -318,7 +315,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& ! find imager sensorindex sensorindex_imager = 0 - iasing_cads = .true. + iasing_cads = .true. !JAJ to be removed if ( iasing_cads .and. imager_coeff ) then if ( sc(2)%sensor_id(1:7) == 'metimag' ) then sensorindex_imager = 2 @@ -474,7 +471,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& ! Check field of view (FOVN) and satellite zenith angle (SAZA) iscn = nint(linele(2)) ! scan line - if( ifov <= 0 .or. ifov > 16) then + if( ifov < 1 .or. ifov > 16) then write(6,*)'READ_IASI-NG: ### ERROR IN READING ', senname, ' BUFR DATA:', & ' STRANGE OBS INFO(FOVN,SLNM):', ifov, iscn cycle read_loop @@ -565,7 +562,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& endif ! Increment nread counter by satinfo_nchan - nread = nread + satinfo_nchan +!JAJ nread = nread + satinfo_nchan crit0 = 0.01_r_kind if( llll > 1 ) crit0 = crit0 + r100 * real(llll,r_kind) @@ -783,6 +780,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& call finalcheck(one,crit1,itx,iuse) endif if(.not. iuse)cycle read_loop + nread = nread + 1 ! Read the imager cluster information for the Cloud and Aerosol Detection Software. ! Only channels 4 and 5 are used. @@ -790,10 +788,9 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& if ( iasing_cads ) then call ufbseq(lnbufr,imager_info,123,7,iret,'IASICSSQ') 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(:)) @@ -813,49 +810,63 @@ subroutine read_iasing(mype,val_iasing,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 18 and 19 are zero, do not compute statistics for the cluster + if ( imager_cluster_size(i) > zero .and. imager_info(105,i) > zero .and. imager_info(111,i) > zero ) then + data_all(maxinfo+j,itx) = imager_cluster_size(i) ! Imager cluster fraction - iexponent = -(nint(imager_info(104,i)) -11 ) ! channel 4 radiance for each cluster. - imager_info(105,i) = imager_info(105,i) * imager_conversion(1) * (ten ** iexponent) + iexponent = -(nint(imager_info(104,i)) -11 ) ! channel 4 radiance for each cluster. + imager_info(105,i) = imager_info(105,i) * imager_conversion(1) * (ten ** iexponent) - iexponent = -(nint(imager_info(106,i)) -11 ) ! channel 4 radiance std dev for each cluster. - imager_info(107,i) = imager_info(107,i) * imager_conversion(1) * (ten ** iexponent) + iexponent = -(nint(imager_info(106,i)) -11 ) ! channel 4 radiance std dev for each cluster. + imager_info(107,i) = imager_info(107,i) * imager_conversion(1) * (ten ** iexponent) - call crtm_planck_temperature(sensorindex_imager,18,imager_info(105,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,18,imager_info(105,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(110,i)) -11 ) ! channel 5 radiance for each cluster - imager_info(111,i) = imager_info(111,i) * imager_conversion(2) * (ten ** iexponent) + iexponent = -(nint(imager_info(110,i)) -11 ) ! channel 5 radiance for each cluster + imager_info(111,i) = imager_info(111,i) * imager_conversion(2) * (ten ** iexponent) - iexponent = -(nint(imager_info(112,i)) -11 ) ! channel 5 radiance std dev for each cluser. - imager_info(113,i) = imager_info(113,i) * imager_conversion(2) * (ten ** iexponent) + iexponent = -(nint(imager_info(112,i)) -11 ) ! channel 5 radiance std dev for each cluser. + imager_info(113,i) = imager_info(113,i) * imager_conversion(2) * (ten ** iexponent) - call crtm_planck_temperature(sensorindex_imager,19,imager_info(110,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,19,imager_info(110,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(105,:)) ! Channel 4 radiance cluster average + imager_mean(1) = sum(imager_cluster_size(:) * imager_info(105,:)) ! Channel 18 radiance cluster average imager_std_dev(1) = sum(imager_cluster_size(:) * (imager_info(105,:)**2 + imager_info(107,:)**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(111,:)) ! Channel 5 radiance cluster average + imager_std_dev(1) = sqrt(max(imager_std_dev(1),zero)) ! Channel 18 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 18 average BT + imager_std_dev(1) = imager_std_dev(1) - imager_mean(1) ! Channel 18 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(111,:)) ! Channel 19 radiance cluster average imager_std_dev(2) = sum(imager_cluster_size(:) * (imager_info(111,:)**2 + imager_info(113,:)**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 19 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 19 average BT + imager_std_dev(2) = imager_std_dev(2) - imager_mean(2) ! Channel 19 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 ! iasing_cads = .true. ! From 714c0a6c042db1fcff93b2e5c310a78e626dcfba Mon Sep 17 00:00:00 2001 From: wx20jjung Date: Mon, 23 Sep 2024 12:12:40 +0000 Subject: [PATCH 06/16] Some changes to read_iasing.f90 --- src/gsi/read_iasing.f90 | 23 ++++++++++------------- 1 file changed, 10 insertions(+), 13 deletions(-) diff --git a/src/gsi/read_iasing.f90 b/src/gsi/read_iasing.f90 index ad3611de21..b0eb521c7d 100644 --- a/src/gsi/read_iasing.f90 +++ b/src/gsi/read_iasing.f90 @@ -562,7 +562,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& endif ! Increment nread counter by satinfo_nchan -!JAJ nread = nread + satinfo_nchan + nread = nread + satinfo_nchan crit0 = 0.01_r_kind if( llll > 1 ) crit0 = crit0 + r100 * real(llll,r_kind) @@ -726,22 +726,21 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& !JAJ !$omp parallel do schedule(dynamic,1) private(i,sc_chan,bufr_chan,radiance) channel_loop: do i=1,satinfo_nchan + sc_chan = sc_index(i) + if (bufr_index(i) == 0 ) cycle channel_loop bufr_chan = bufr_index(i) - if (bufr_chan > 0 ) then ! check that channel number is within reason - if (( allchan(2,bufr_chan) > zero .and. allchan(2,bufr_chan) < 99999._r_kind)) then ! radiance bounds - radiance = allchan(2,bufr_chan)*scalef(bufr_chan) - sc_chan = sc_index(i) - call crtm_planck_temperature(sensorindex_iasing,sc_chan,radiance,temperature(bufr_chan)) + if (( allchan(2,bufr_chan) > zero .and. allchan(2,bufr_chan) < 99999._r_kind)) then ! radiance bounds + radiance = allchan(2,bufr_chan)*scalef(bufr_chan) + call crtm_planck_temperature(sensorindex_iasing,sc_chan,radiance,temperature(bufr_chan)) !JAJ if (( temperature(bufr_chan) > 220.0_r_kind .and. temperature(bufr_chan) < 300.0_r_kind) .and. & !JAJ (bufr_chan > 500 .and. bufr_chan < 5000)) then !JAJ write(*,'(a10,1x,i6,1x,e15.3,1x,f7.2,1x,f10.5,1x,f10.5)') 'JAJ value ',bufr_chan,allchan(2,bufr_chan),temperature(bufr_chan), & !JAJ dlat_earth_deg, dlon_earth_deg !JAJ endif - else - temperature(bufr_chan) = tbmin - endif + else + temperature(bufr_chan) = tbmin end if end do channel_loop @@ -751,7 +750,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& if ( bufr_index(i) == 0 ) cycle skip_loop bufr_chan = bufr_index(i) if(temperature(bufr_chan) <= tbmin .or. temperature(bufr_chan) > tbmax ) then - temperature(bufr_chan) = min(tbmax,max(tbmin,temperature(bufr_chan))) + temperature(bufr_chan) = tbmin !write(*,'(a10,1x,i6,1x,e7.3,1x,f9.5,1x,f9.5)') 'JAJ bad value ',bufr_chan,allchan(2,bufr_chan),scalef(bufr_chan) ! write(*,'(a14,1x,i6,1x,e12.3,1x,f10.5,1x,f10.5)') 'JAJ bad value ',bufr_chan,allchan(2,bufr_chan),dlat_earth_deg, dlon_earth_deg if(iuse_rad(ioff+i) >= 0)iskip = iskip + 1 @@ -780,7 +779,6 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& call finalcheck(one,crit1,itx,iuse) endif if(.not. iuse)cycle read_loop - nread = nread + 1 ! Read the imager cluster information for the Cloud and Aerosol Detection Software. ! Only channels 4 and 5 are used. @@ -931,9 +929,8 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& ! Put satinfo defined channel temperatures into data array do l=1,satinfo_nchan ! Prevent out of bounds reference from temperature - if ( bufr_index(l) == 0 ) cycle i = bufr_index(l) - if(i /= 0)then + if(bufr_index(l) /= 0) then data_all(l+nreal,itx) = temperature(i) ! brightness temerature else data_all(l+nreal,itx) = tbmin From cd3d7e859035337bd136e760ff6518fea10989ec Mon Sep 17 00:00:00 2001 From: wx20jjung Date: Mon, 30 Sep 2024 11:33:43 +0000 Subject: [PATCH 07/16] Some code cleanup in read_iasing.f90 --- src/gsi/read_iasing.f90 | 14 +++----------- 1 file changed, 3 insertions(+), 11 deletions(-) diff --git a/src/gsi/read_iasing.f90 b/src/gsi/read_iasing.f90 index b0eb521c7d..eec4868ce2 100644 --- a/src/gsi/read_iasing.f90 +++ b/src/gsi/read_iasing.f90 @@ -315,7 +315,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& ! find imager sensorindex sensorindex_imager = 0 - iasing_cads = .true. !JAJ to be removed +!JAJ iasing_cads = .true. !JAJ to be removed if ( iasing_cads .and. imager_coeff ) then if ( sc(2)%sensor_id(1:7) == 'metimag' ) then sensorindex_imager = 2 @@ -681,7 +681,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& endif end do -!JAJ NEED TO CHECK THIS ONB +!JAJ NEED TO CHECK THIS ONE ! Read IASI-NG channel number(CHNM) and radiance (SCRA) call ufbseq(lnbufr,allchan,2,bufr_nchan,iret,'I1CRSQ') jstart=1 @@ -724,7 +724,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& cycle read_loop endif -!JAJ !$omp parallel do schedule(dynamic,1) private(i,sc_chan,bufr_chan,radiance) +!$omp parallel do schedule(dynamic,1) private(i,sc_chan,bufr_chan,radiance) channel_loop: do i=1,satinfo_nchan sc_chan = sc_index(i) if (bufr_index(i) == 0 ) cycle channel_loop @@ -733,12 +733,6 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& if (( allchan(2,bufr_chan) > zero .and. allchan(2,bufr_chan) < 99999._r_kind)) then ! radiance bounds radiance = allchan(2,bufr_chan)*scalef(bufr_chan) call crtm_planck_temperature(sensorindex_iasing,sc_chan,radiance,temperature(bufr_chan)) -!JAJ if (( temperature(bufr_chan) > 220.0_r_kind .and. temperature(bufr_chan) < 300.0_r_kind) .and. & -!JAJ (bufr_chan > 500 .and. bufr_chan < 5000)) then -!JAJ write(*,'(a10,1x,i6,1x,e15.3,1x,f7.2,1x,f10.5,1x,f10.5)') 'JAJ value ',bufr_chan,allchan(2,bufr_chan),temperature(bufr_chan), & -!JAJ dlat_earth_deg, dlon_earth_deg -!JAJ endif - else temperature(bufr_chan) = tbmin end if @@ -751,8 +745,6 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& bufr_chan = bufr_index(i) if(temperature(bufr_chan) <= tbmin .or. temperature(bufr_chan) > tbmax ) then temperature(bufr_chan) = tbmin - !write(*,'(a10,1x,i6,1x,e7.3,1x,f9.5,1x,f9.5)') 'JAJ bad value ',bufr_chan,allchan(2,bufr_chan),scalef(bufr_chan) -! write(*,'(a14,1x,i6,1x,e12.3,1x,f10.5,1x,f10.5)') 'JAJ bad value ',bufr_chan,allchan(2,bufr_chan),dlat_earth_deg, dlon_earth_deg if(iuse_rad(ioff+i) >= 0)iskip = iskip + 1 endif end do skip_loop From db33a1ba051f8765742858c9307d0fcd9dc67e8e Mon Sep 17 00:00:00 2001 From: "jim.jung" Date: Fri, 18 Oct 2024 13:44:08 +0000 Subject: [PATCH 08/16] Changes to read_iasing.f90 --- src/gsi/read_iasing.f90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/gsi/read_iasing.f90 b/src/gsi/read_iasing.f90 index b0eb521c7d..6141795629 100644 --- a/src/gsi/read_iasing.f90 +++ b/src/gsi/read_iasing.f90 @@ -185,7 +185,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& ! Set standard parameters character(8),parameter:: fov_flag="crosstrk" - integer(i_kind),parameter:: sfc_channel=2539 + integer(i_kind),parameter:: sfc_channel=2541 integer(i_kind),parameter:: ichan=-999 ! fov-based surface code is not channel specific for iasing real(r_kind),parameter:: expansion=one ! exansion factor for fov-based surface code. ! use one for ir sensors. @@ -681,7 +681,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& endif end do -!JAJ NEED TO CHECK THIS ONB +!JAJ NEED TO CHECK THIS ONE ! Read IASI-NG channel number(CHNM) and radiance (SCRA) call ufbseq(lnbufr,allchan,2,bufr_nchan,iret,'I1CRSQ') jstart=1 @@ -724,7 +724,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& cycle read_loop endif -!JAJ !$omp parallel do schedule(dynamic,1) private(i,sc_chan,bufr_chan,radiance) +!$omp parallel do schedule(dynamic,1) private(i,sc_chan,bufr_chan,radiance) channel_loop: do i=1,satinfo_nchan sc_chan = sc_index(i) if (bufr_index(i) == 0 ) cycle channel_loop @@ -759,7 +759,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& if(iskip > 0)then if(print_verbose)write(6,*) ' READ_IASI-NG: iskip > 0 ',iskip - cycle read_loop +!JAJ add cycle read_loop end if ! crit1=crit1 + ten*real(iskip,r_kind) From 9446d500f7e683133b8753d6aa2ac25d21e407f9 Mon Sep 17 00:00:00 2001 From: "jim.jung" Date: Sat, 2 Nov 2024 19:51:35 +0000 Subject: [PATCH 09/16] Modified read_iasing to support the use_edges flag and remove the scan angle calculation. --- src/gsi/read_iasing.f90 | 234 +++++++++++++++++----------------------- 1 file changed, 98 insertions(+), 136 deletions(-) diff --git a/src/gsi/read_iasing.f90 b/src/gsi/read_iasing.f90 index be7bfd3cd1..d857206f40 100644 --- a/src/gsi/read_iasing.f90 +++ b/src/gsi/read_iasing.f90 @@ -62,7 +62,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& use satthin, only: radthin_time_info,tdiff2crit use obsmod, only: time_window_max use radinfo, only:iuse_rad,nuchan,nusis,jpch_rad,crtm_coeffs_path,use_edges, & - radedge1,radedge2,radstart,radstep + radedge1,radedge2 use crtm_module, only: success, & crtm_kind => fp use crtm_planck_functions, only: crtm_planck_temperature @@ -117,10 +117,9 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& real(r_double),dimension(13) :: allspot real(r_double),allocatable,dimension(:,:) :: allchan real(r_double),dimension(3,4):: cscale - real(r_double),dimension(7):: cloud_frac + real(r_double),dimension(7):: cluster_frac integer(i_kind) :: bufr_size - real(r_kind) :: step, start,step_adjust character(len=8) :: subset character(len=4) :: senname character(len=80) :: allspotlist @@ -140,7 +139,6 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& real(r_kind) :: piece real(r_kind) :: rsat, dlon, dlat real(r_kind) :: dlon_earth,dlat_earth,dlon_earth_deg,dlat_earth_deg - real(r_kind) :: lza, lzaest,sat_height_ratio real(r_kind) :: pred, crit1, dist1 real(r_kind) :: sat_zenang real(crtm_kind) :: radiance @@ -156,21 +154,21 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& real(r_kind) cdist,disterr,disterrmax,dlon00,dlat00 logical :: outside,iuse,assim,valid - logical :: iasing,quiet,cloud_info + logical :: quiet,cloud_info integer(i_kind) :: ifov, ifor, istep, ipos, instr, iscn, ioff, sensorindex_iasing integer(i_kind) :: i, j, l, iskip, ifovn, bad_line, ksatid, kidsat, llll integer(i_kind) :: nreal, isflg integer(i_kind) :: itx, k, nele, itt, n - integer(i_kind):: iexponent,maxinfo, bufr_nchan, dval_info - integer(i_kind):: idomsfc(1) - integer(i_kind):: ntest - integer(i_kind):: error_status, irecx,ierr - integer(i_kind):: radedge_min, radedge_max - integer(i_kind) :: subset_start, subset_end, satinfo_nchan, sc_chan, bufr_chan - integer(i_kind) :: sfc_channel_index - integer(i_kind),allocatable, dimension(:) :: channel_number, sc_index, bufr_index - integer(i_kind),allocatable, dimension(:) :: bufr_chan_test + integer(i_kind) :: iexponent,maxinfo, bufr_nchan, dval_info + integer(i_kind) :: idomsfc(1) + integer(i_kind) :: ntest + integer(i_kind) :: error_status, irecx,ierr + integer(i_kind) :: radedge_min, radedge_max + integer(i_kind) :: subset_start, subset_end, satinfo_nchan, sc_chan, bufr_chan + integer(i_kind) :: sfc_channel_index + integer(i_kind),allocatable, dimension(:) :: channel_number, sc_index, bufr_index + integer(i_kind),allocatable, dimension(:) :: bufr_chan_test character(len=20),allocatable, dimension(:):: sensorlist ! Imager cluster information for CADS @@ -187,13 +185,12 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& character(8),parameter:: fov_flag="crosstrk" integer(i_kind),parameter:: sfc_channel=2541 integer(i_kind),parameter:: ichan=-999 ! fov-based surface code is not channel specific for iasing - real(r_kind),parameter:: expansion=one ! exansion factor for fov-based surface code. + real(r_kind),parameter:: expansion=one ! expansion factor for fov-based surface code. ! use one for ir sensors. real(r_kind),parameter:: R90 = 90._r_kind real(r_kind),parameter:: R360 = 360._r_kind - real(r_kind),parameter:: tbmin = 50._r_kind + real(r_kind),parameter:: tbmin = 50._r_kind real(r_kind),parameter:: tbmax = 550._r_kind - real(r_kind),parameter:: earth_radius = 6371000._r_kind integer(i_kind),parameter :: ilon = 3 integer(i_kind),parameter :: ilat = 4 real(r_kind) :: ptime,timeinflat,crit0 @@ -215,7 +212,6 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& ndata = 0 nodata = 0 - iasing = obstype == 'iasi-ng' bad_line=-1 @@ -223,7 +219,12 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& call gsi_nstcoupler_skindepth(obstype, zob) ! get penetration depth (zob) for the obstype endif - if (jsatid == 'metop-sg-a1') kidsat=24 + if (jsatid == 'metop-sg-a1') then + kidsat=24 + else + write(*,*) 'READ_IASI-NG: Unrecognized value for jsatid '//jsatid//': RETURNING' + return + endif ! write(6,*)'READ_IASI-NG: mype, mype_root,mype_sub, npe_sub,mpi_comm_sub', & ! mype, mype_root,mype_sub,mpi_comm_sub @@ -241,8 +242,6 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& if (trim(nusis(i))==trim(sis)) then ioff = min(ioff,i) ! iasi-ng offset if (subset_start == 0) then -!JAJ step = radstep(i) -!JAJ start = radstart(i) if (radedge1(i)/=-1 .and. radedge2(i)/=-1) then radedge_min=radedge1(i) radedge_max=radedge2(i) @@ -259,14 +258,13 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& allocate(bufr_index(satinfo_nchan)) ioff = ioff - 1 - step_adjust = 0.625_r_kind ! If all channels of a given sensor are set to monitor or not ! assimilate mode (iuse_rad<1), reset relative weight to zero. ! We do not want such observations affecting the relative ! weighting between observations within a given thinning group. if (.not.assim) val_iasing=zero - if (mype_sub==mype_root)write(6,*)'READ_IASI-NG: iasi-ng offset ',ioff + if (mype_sub==mype_root) write(6,*)'READ_IASI-NG: iasi-ng offset ',ioff senname = 'IASI-NG' @@ -304,8 +302,8 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& ! find IASI-NG sensorindex sensorindex_iasing = 0 + if ( sc(1)%sensor_id(1:7) == 'iasi-ng' .or. sc(1)%sensor_id == '999') then !JAJ if ( sc(1)%sensor_id(1:7) == 'iasi-ng' ) then - if ( sc(1)%sensor_id == '999' ) then sensorindex_iasing = 1 else write(6,*)'READ_IASI-NG: ***ERROR*** sensorindex_iasi-ng not set NO IASI-NG DATA USED' @@ -315,7 +313,6 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& ! find imager sensorindex sensorindex_imager = 0 -!JAJ iasing_cads = .true. !JAJ to be removed if ( iasing_cads .and. imager_coeff ) then if ( sc(2)%sensor_id(1:7) == 'metimag' ) then sensorindex_imager = 2 @@ -345,7 +342,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& ! Calculate parameters needed for FOV-based surface calculation. if (isfcalc==1)then - instr=18 + instr=18 !This is actually IASI. Must be changed for IASI-NG if used. call instrument_init(instr, jsatid, expansion, valid) if (.not. valid) then if (assim) then @@ -406,11 +403,12 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& infile2=trim(infile)//'_db' ! Set bufr subset names based on type of data to read end if -! Open BUFR file +! Open BUFR file open(lnbufr,file=trim(infile2),form='unformatted',status='old',iostat=ierr) if(ierr /= 0) cycle ears_db_loop -! Open BUFR table + +! Open BUFR table call openbf(lnbufr,'IN',lnbufr) call datelen(10) @@ -425,7 +423,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& read_loop: do while (ireadsb(lnbufr)==0) -! Get the size of the channels and radiance (allchan) array +! Get the size of the channels and radiance (allchan) array call ufbint(lnbufr,crchn_reps,1,1,iret,'(I1CRSQ)') bufr_nchan = int(crchn_reps) @@ -440,16 +438,16 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& bufr_chan_test(:)=0 endif ! allocation if -! Read IASI-NG FOV information +! Read IASI-NG FOV information call ufbint(lnbufr,linele,5,1,iret,'FOVN SLNM INGGQF SELV SAID') -! Extract satellite id. If not the one we want, read next subset +! Extract satellite id. If not the one we want, read next subset ksatid=nint(linele(5)) if(ksatid /= kidsat) cycle read_loop if ( linele(3) /= zero) cycle read_loop ! problem with profile (INGGQF) -! zenith angle/scan spot mismatch, reject entire line +! Zenith angle/scan spot mismatch, reject entire line if ( bad_line == nint(linele(2))) then cycle read_loop else @@ -457,30 +455,16 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& endif ifov = nint(linele(1)) ! field of view - -! IASI-NG fov ranges from 1 to 112. Current angle dependent bias -! correction has a maximum of 90 scan positions. Geometry -! of IASI-NG scan allows us to remap 1-112 to 1-56. Variable -! ifovn below contains the remapped IASI-NG fov. This value is -! passed on to and used in setuprad -!JAJ ifovn = (ifov-1)/2 + 1 - -! Remove data on edges -!JAJ if (.not. use_edges .and. & -!JAJ (ifovn < radedge_min .OR. ifovn > radedge_max )) cycle read_loop - -! Check field of view (FOVN) and satellite zenith angle (SAZA) - iscn = nint(linele(2)) ! scan line - if( ifov < 1 .or. ifov > 16) then + if ( ifov < 1 .or. ifov > 224 ) then ! field of view out of range write(6,*)'READ_IASI-NG: ### ERROR IN READING ', senname, ' BUFR DATA:', & - ' STRANGE OBS INFO(FOVN,SLNM):', ifov, iscn + ' BAD FIELD OF VIEW:', ifov cycle read_loop endif call ufbint(lnbufr,allspot,13,1,iret,allspotlist) if(iret /= 1) cycle read_loop -! Check observing position +! Check observing position dlat_earth = allspot(8) ! latitude dlon_earth = allspot(9) ! longitude if( abs(dlat_earth) > R90 .or. abs(dlon_earth) > R360 .or. & @@ -490,7 +474,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& cycle read_loop endif -! Retrieve observing position +! Retrieve observing position if(dlon_earth >= R360)then dlon_earth = dlon_earth - R360 else if(dlon_earth < ZERO)then @@ -502,7 +486,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& dlat_earth = dlat_earth * deg2rad dlon_earth = dlon_earth * deg2rad -! If regional, map obs lat,lon to rotated grid. +! If regional, map obs lat,lon to rotated grid. if(regional)then ! Convert to rotated coordinate. dlon centered on 180 (pi), @@ -517,12 +501,11 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& disterr=acos(cdist)*rad2deg disterrmax=max(disterrmax,disterr) end if - -! Check to see if in domain. outside=.true. if dlon_earth, -! dlat_earth outside domain, =.false. if inside + +! Check to see if in domain. outside=.true. if dlon_earth and dlat_earth outside domain =.false. if inside if(outside) cycle read_loop -! Global case +! Global case else dlat = dlat_earth dlon = dlon_earth @@ -530,7 +513,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& call grdcrd1(dlon,rlons,nlon,1) endif -! Check obs time +! Check obs time idate5(1) = nint(allspot(2)) ! year idate5(2) = nint(allspot(3)) ! month idate5(3) = nint(allspot(4)) ! day @@ -549,7 +532,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& endif -! Retrieve obs time +! Retrieve obs time call w3fs21(idate5,nmind) t4dv = (real(nmind-iwinbgn,r_kind) + real(allspot(7),r_kind)*r60inv)*r60inv ! add in seconds sstime = real(nmind,r_kind) + real(allspot(7),r_kind)*r60inv ! add in seconds @@ -561,7 +544,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& if (abs(tdiff)>twind) cycle read_loop endif -! Increment nread counter by satinfo_nchan +! Increment nread counter by satinfo_nchan nread = nread + satinfo_nchan crit0 = 0.01_r_kind @@ -572,52 +555,38 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& if(.not. iuse)cycle read_loop -! Observational info +! Observational info sat_zenang = allspot(10) ! satellite zenith angle -! Check satellite zenith angle (SAZA) +! Check satellite zenith angle (SAZA) if(sat_zenang > 90._r_kind ) then write(6,*)'READ_IASI-NG: ### ERROR IN READING ', senname, ' BUFR DATA:', & ' STRANGE OBS INFO(FOVN,SLNM,SAZA,BEARAZ):', ifov, iscn, allspot(10),allspot(11) cycle read_loop endif -! Remove this code when the field-of-view spans 1-224 JAJ. - ifov = ( 59 - nint(sat_zenang)) / 2 - ifov = min( 28, max(1,ifov)) - -! There are 16 fields-of-view within each field-of-reagard, a 4X4 matrix. -! | 1 2 3 4 | -! | 5 6 7 8 | -! | 9 10 11 12 | -! | 13 14 15 16 | - -! There are 14 fields-of-regard in each scan line. So, there are 56 unique positions in the scan line. -! To determine the scan position: - ifor = ((ifov-1) / 16) + 1 ! Determine field-of-regard - istep = ifov - (ifor * 16) ! Determine field-of-view within field-of-regard - ipos = (ifor * 4) + mod(istep,4) + 1 ! Determine position of field-of-view within scan line + +! There are 16 fields-of-view within each field-of-reagard, a 4X4 matrix. +! | 1 2 3 4 | +! | 5 6 7 8 | +! | 9 10 11 12 | +! | 13 14 15 16 | + +! Add this code when field-of-view spans 1-224 JAJ. +! There are 14 fields-of-regard in each scan line. So, there are 56 unique positions in the scan line. +! To determine the scan position: + ifor = (ifov-1) / 16 ! Determine field-of-regard + istep = ifov - ((ifor) * 16) ! Determine field-of-view within field-of-regard + ipos = (ifor * 4) + mod(istep,4) ! Determine position of field-of-view within scan line ! Remove data on edges - if (.not. use_edges .and. & - (ipos < radedge_min .or. ipos > radedge_max)) cycle read_loop - - if ( ipos <= 28 ) sat_zenang = - sat_zenang - -!JAJ to be removed -! Estimate scan angle and compare to zenith angle -!JAJ lza = ((start + real(ipos)* step) * deg2rad -!JAJ sat_height_ratio = (earth_radius + linele(4)) / earth_radius -!JAJ lzaest = asin(sat_height_ratio * sin(lza)) * rad2deg -!JAJ if (abs(sat_zenang - lzaest) > one) then -!JAJ write(6,*)' READ_IASI-NG WARNING uncertainty in lza ', & -!JAJ lza*rad2deg,sat_zenang,sis,ifov,start,step,allspot(11),allspot(12),allspot(13) -!JAJ bad_line = iscn -!JAJ cycle read_loop -!JAJ endif - -! "Score" observation. We use this information to identify "best" obs -! Locate the observation on the analysis grid. Get sst and land/sea/ice -! mask. + if (.not. use_edges .and. & + (ipos < radedge_min .or. ipos > radedge_max)) cycle read_loop + +! Set sign of satellite zenith angle + if (ipos <= 28) sat_zenang = -sat_zenang + +! "Score" observation. We use this information to identify "best" obs +! Locate the observation on the analysis grid. Get sst and land/sea/ice mask. ! isflg - surface flag ! 0 sea ! 1 land @@ -625,13 +594,13 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& ! 3 snow ! 4 mixed -! When using FOV-based surface code, must screen out obs with bad fov numbers. +! When using FOV-based surface code, must screen out obs with bad fov numbers. if (isfcalc == 1) then call fov_check(ifov,instr,ichan,valid) if (.not. valid) cycle read_loop -! When isfcalc is set to one, calculate surface fields using size/shape of fov. -! Otherwise, use bilinear interpolation. +! When isfcalc is set to one, calculate surface fields using size/shape of fov. +! Otherwise, use bilinear interpolation. call deter_sfc_fov(fov_flag,ifov,instr,ichan,real(allspot(11),r_kind),dlat_earth_deg, & dlon_earth_deg,expansion,t4dv,isflg,idomsfc(1), & @@ -641,19 +610,19 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& ts,tsavg,vty,vfr,sty,stp,sm,sn,zz,ff10,sfcr) endif -! Set common predictor parameters +! Set common predictor parameters crit1 = crit1 + rlndsea(isflg) call checkob(one,crit1,itx,iuse) if(.not. iuse)cycle read_loop -! Clear Amount (percent clear) -! Compute "score" for observation. All scores>=0.0. Lowest score is "best" +! Compute "score" for observation. All scores>=0.0. Lowest score is "best" pred = r100 cloud_info = .false. - call ufbrep(lnbufr,cloud_frac,1,7,iret,'ASCS') - if (iret == 7 .and. cloud_frac(1) <= r100 .and. cloud_frac(1) >= zero) then - pred = r100 - cloud_frac(1) + call ufbrep(lnbufr,cluster_frac,1,7,iret,'ASCS') +! Use size of warmest (first) cluster as thinning criteria. + if (iret == 7 .and. cluster_frac(1) <= r100 .and. cluster_frac(1) >= zero) then + pred = r100 - cluster_frac(1) cloud_info = .true. endif @@ -667,11 +636,11 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& cycle read_loop end if -! The scaling factors are as follows, cscale(1) is the start channel number, +! The scaling factors are as follows, cscale(1) is the start channel number, ! cscale(2) is the end channel number, ! cscale(3) is the exponent scaling factor -! In our case, there are 4 groups of cscale (dimension :: cscale(3,4)) -! The units are W/m2..... you need to convert to mW/m2.... (subtract 5 from cscale(3,:) +! In our case, there are 4 groups of cscale (dimension :: cscale(3,4)) +! The units are W/m2..... you need to convert to mW/m2.... (subtract 5 from cscale(3,:) do i=1,4 ! convert exponent scale factor to int and change units if(cscale(3,i) < bmiss) then iexponent = -(nint(cscale(3,i)) - 5) @@ -681,8 +650,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& endif end do -!JAJ NEED TO CHECK THIS ONE -! Read IASI-NG channel number(CHNM) and radiance (SCRA) +! Read IASI-NG channel number(CHNM) and radiance (SCRA) call ufbseq(lnbufr,allchan,2,bufr_nchan,iret,'I1CRSQ') jstart=1 scalef=one @@ -702,8 +670,8 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& cycle read_loop endif -! Coordinate bufr channels with satinfo file channels -! If this is the first time or a change in the bufr channels is detected, sync with satinfo file +! Coordinate bufr channels with satinfo file channels +! If this is the first time or a change in the bufr channels is detected, sync with satinfo file if (ANY(int(allchan(1,:)) /= bufr_chan_test(:))) then sfc_channel_index = 0 bufr_index(:) = 0 @@ -738,7 +706,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& end if end do channel_loop -! Check for reasonable temperature values +! Check for reasonable temperature values iskip = 0 skip_loop: do i=1,satinfo_nchan if ( bufr_index(i) == 0 ) cycle skip_loop @@ -750,21 +718,21 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& end do skip_loop if(iskip > 0)then - if(print_verbose)write(6,*) ' READ_IASI-NG: iskip > 0 ',iskip -!JAJ add cycle read_loop + if(print_verbose)write(6,*) ' READ_IASI-NG: iskip > 0 ',iskip + cycle read_loop end if ! crit1=crit1 + ten*real(iskip,r_kind) -! If the surface channel exists (~960.0 cm-1) and the imager cloud information is missing, use an -! estimate of the surface temperature to determine if the profile may be clear. +! If the surface channel exists (~960.0 cm-1) and the imager cloud information is missing, use an +! estimate of the surface temperature to determine if the profile may be clear. if (.not. cloud_info) then pred = tsavg*0.98_r_kind - temperature(sfc_channel_index) pred = max(pred,zero) crit1=crit1 + pred endif -! Map obs to grids +! Map obs to grids if (pred == zero) then call finalcheck(dist1,crit1,itx,iuse) else @@ -772,8 +740,8 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& endif if(.not. iuse)cycle read_loop -! Read the imager cluster information for the Cloud and Aerosol Detection Software. -! Only channels 4 and 5 are used. +! Read the imager cluster information for the Cloud and Aerosol Detection Software. +! Only channels 18 and 19 are used. if ( iasing_cads ) then call ufbseq(lnbufr,imager_info,123,7,iret,'IASICSSQ') @@ -787,20 +755,19 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& imager_conversion(1) = one / (sc(sensorindex_imager)%wavenumber(18) **2) imager_conversion(2) = one / (sc(sensorindex_imager)%wavenumber(19) **2) -! Order clusters from largest (1) to smallest (7) +! Order clusters from largest (1) to smallest (7) imager_cluster_sort: do i=1,7 j = maxloc(imager_cluster_size,dim=1,mask=imager_cluster_flag) imager_cluster_index(i) = j imager_cluster_flag(j) = .FALSE. end do imager_cluster_sort -! Convert from radiance to brightness temperature for mean and standard deviation used by CADS. -! Imager cluster info added to data_all array - +! Convert from radiance to brightness temperature for mean and standard deviation used by CADS. +! Imager cluster info added to data_all array imager_cluster_info: do j=1,7 i = imager_cluster_index(j) -! If the cluster size, or radiance values of channel 18 and 19 are zero, do not compute statistics for the cluster +! If the cluster size, or radiance values of channel 18 and 19 are zero, do not compute statistics for the cluster if ( imager_cluster_size(i) > zero .and. imager_info(105,i) > zero .and. imager_info(111,i) > zero ) then data_all(maxinfo+j,itx) = imager_cluster_size(i) ! Imager cluster fraction @@ -830,7 +797,6 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& end do imager_cluster_info ! Compute cluster averages for each channel - imager_mean(1) = sum(imager_cluster_size(:) * imager_info(105,:)) ! Channel 18 radiance cluster average imager_std_dev(1) = sum(imager_cluster_size(:) * (imager_info(105,:)**2 + imager_info(107,:)**2)) - imager_mean(1)**2 imager_std_dev(1) = sqrt(max(imager_std_dev(1),zero)) ! Channel 18 radiance RMSE @@ -859,9 +825,8 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& data_all(maxinfo+1 : maxinfo+cads_info,itx) = zero endif endif ! iasing_cads = .true. -! -! interpolate NSST variables to Obs. location and get dtw, dtc, tz_tr -! + +! Interpolate NSST variables to Obs. location and get dtw, dtc, tz_tr if ( nst_gsi > 0 ) then tref = ts(0) dtw = zero @@ -879,8 +844,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& data_all(4,itx) = dlat ! grid relative latitude data_all(5,itx) = sat_zenang*deg2rad ! satellite zenith angle (rad) data_all(6,itx) = allspot(11) ! satellite azimuth angle (deg) -!JAJ data_all(7,itx) = lza ! look angle (rad) - data_all(7,itx) = sat_zenang*deg2rad ! look angle (rad) + data_all(7,itx) = sat_zenang*deg2rad ! look angle (rad) data_all(8,itx) = ifov ! fov number data_all(9,itx) = allspot(12) ! solar zenith angle (deg) data_all(10,itx)= allspot(13) ! solar azimuth angle (deg) @@ -918,7 +882,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& data_all(maxinfo+cads_info+dval_info+4,itx) = tz_tr ! d(Tz)/d(Tr) endif -! Put satinfo defined channel temperatures into data array +! Put satinfo defined channel temperatures into data array do l=1,satinfo_nchan ! Prevent out of bounds reference from temperature i = bufr_index(l) @@ -949,7 +913,6 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& ! If multiple tasks read input bufr file, allow each tasks to write out ! information it retained and then let single task merge files together - call combine_radobs(mype_sub,mype_root,npe_sub,mpi_comm_sub,& nele,itxmax,nread,ndata,data_all,score_crit,nrec) @@ -957,9 +920,8 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& ! and write out data to scratch file for further processing. if (mype_sub==mype_root.and.ndata>0) then -! Identify "bad" observation (unreasonable brightness temperatures). -! Update superobs sum according to observation location - +! Identify "bad" observation (unreasonable brightness temperatures). +! Update superobs sum according to observation location do n=1,ndata do i=1,satinfo_nchan if(data_all(i+nreal,n) > tbmin .and. & @@ -974,7 +936,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& end do end if -! Write final set of "best" observations to output file +! Write final set of "best" observations to output file call count_obs(ndata,nele,ilat,ilon,data_all,nobs) write(lunout) obstype,sis,nreal,satinfo_nchan,ilat,ilon write(lunout) ((data_all(k,n),k=1,nele),n=1,ndata) From f5027badcf07e801719d91e69a457ea4d98ef24b Mon Sep 17 00:00:00 2001 From: James Jung Date: Tue, 5 Nov 2024 09:41:06 -0500 Subject: [PATCH 10/16] Update src/gsi/gsimod.F90 Co-authored-by: David Huber <69919478+DavidHuber-NOAA@users.noreply.github.com> --- src/gsi/gsimod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index 99bedda7e9..3c5ed07329 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -1064,7 +1064,7 @@ module gsimod ! ! Flags to use the new IR cloud detection routine. Flag must be set to true to use the new routine. The default ! (no flag or .false.) will use the default. -! airs_cads : use the clod and aerosool detection software for the AIRS instrument +! airs_cads : use the cloud and aerosol detection software for the AIRS instrument ! cris_cads : use the cloud and aerosol detection software for CrIS instruments ! iasi_cads : use the cloud and aerosol detection software for IASI instruments ! iasing_cads: use the cloud and aerosol detection software for IASI instruments From 68037a07221ae4aecf95a2f818a906a59a1f444a Mon Sep 17 00:00:00 2001 From: James Jung Date: Tue, 5 Nov 2024 09:50:38 -0500 Subject: [PATCH 11/16] Update src/gsi/gsimod.F90 Yes, thanks Co-authored-by: David Huber <69919478+DavidHuber-NOAA@users.noreply.github.com> --- src/gsi/gsimod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index 3c5ed07329..2e4a1847ce 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -1067,7 +1067,7 @@ module gsimod ! airs_cads : use the cloud and aerosol detection software for the AIRS instrument ! cris_cads : use the cloud and aerosol detection software for CrIS instruments ! iasi_cads : use the cloud and aerosol detection software for IASI instruments -! iasing_cads: use the cloud and aerosol detection software for IASI instruments +! iasing_cads: use the cloud and aerosol detection software for IASI-NG instruments ! namelist/obsqc/dfact,dfact1,erradar_inflate,tdrerr_inflate,oberrflg,& From e2560566c3973d077a4fbae5dda24bc21f5d8232 Mon Sep 17 00:00:00 2001 From: wx20jjung Date: Tue, 5 Nov 2024 15:24:28 +0000 Subject: [PATCH 12/16] These are changes suggested by Dave Huber. Added comments and removed code that is no longer used. --- src/gsi/crtm_interface.f90 | 4 ++-- src/gsi/gsimod.F90 | 4 ++-- src/gsi/qcmod.f90 | 4 ++-- src/gsi/read_iasing.f90 | 16 ++-------------- 4 files changed, 8 insertions(+), 20 deletions(-) diff --git a/src/gsi/crtm_interface.f90 b/src/gsi/crtm_interface.f90 index 4120a9a3a8..2f2764b827 100644 --- a/src/gsi/crtm_interface.f90 +++ b/src/gsi/crtm_interface.f90 @@ -667,8 +667,8 @@ subroutine init_crtm(init_pass,mype_diaghdr,mype,nchanl,nreal,isis,obstype,radmo error_status = crtm_channelinfo_subset(channelinfo(1), & channel_subset = nuchan(subset_start:subset_end)) -!JAJ the correct one else if (channelinfo(1)%sensor_id(1:7) == 'iasi-ng' .AND. isis(1:7) == 'iasi-ng') then -else if (channelinfo(1)%sensor_id(1:3) == '999' .AND. isis(1:7) == 'iasi-ng') then +!TODO JAJ the correct one else if (channelinfo(1)%sensor_id(1:7) == 'iasi-ng' .AND. isis(1:7) == 'iasi-ng') then +else if (channelinfo(1)%sensor_id(1:3) == '999' .AND. isis(1:7) == 'iasi-ng') then ! TODO To be removed. sensorindex = 1 subset_start = 0 subset_end = 0 diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index 99bedda7e9..8bae3ca613 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -1064,10 +1064,10 @@ module gsimod ! ! Flags to use the new IR cloud detection routine. Flag must be set to true to use the new routine. The default ! (no flag or .false.) will use the default. -! airs_cads : use the clod and aerosool detection software for the AIRS instrument +! airs_cads : use the cloud and aerosol detection software for AIRS instrument ! cris_cads : use the cloud and aerosol detection software for CrIS instruments ! iasi_cads : use the cloud and aerosol detection software for IASI instruments -! iasing_cads: use the cloud and aerosol detection software for IASI instruments +! iasing_cads: use the cloud and aerosol detection software for IASI-NG instruments ! namelist/obsqc/dfact,dfact1,erradar_inflate,tdrerr_inflate,oberrflg,& diff --git a/src/gsi/qcmod.f90 b/src/gsi/qcmod.f90 index f1e81323c8..86207a49fb 100644 --- a/src/gsi/qcmod.f90 +++ b/src/gsi/qcmod.f90 @@ -2112,6 +2112,8 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr,airs ! goessndr - logical flag - if goessndr data - true ! cris - logical flag - if cris data - true ! avhrr - logical flag - if avhrr data - true +! iasi - logical flag - if iasi data - true +! iasing - logical flag - if iasing data - true ! zsges - elevation of guess ! cenlat - latitude of observation ! frac_sea - fraction of grid box covered with water @@ -2434,8 +2436,6 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr,airs ! default compute cloud stats, emc_legacy_cloud_detect else if ( lcloud > 0 ) then -!JAJ if ( lcloud > 0 .and. .not. iasing ) then - do i=1,nchanl ! reject channels with iuse_rad(j)=-1 when they are peaking below the cloud j=ich(i) diff --git a/src/gsi/read_iasing.f90 b/src/gsi/read_iasing.f90 index d857206f40..c20de10524 100644 --- a/src/gsi/read_iasing.f90 +++ b/src/gsi/read_iasing.f90 @@ -43,6 +43,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& ! nrec_start - first subset with useful information ! nrec_start_ears - first ears subset with useful information ! nrec_start_db - first db subset with useful information +! dval_use - logical for using dval ! ! output argument list: ! nread - number of BUFR IASI-NG observations read @@ -50,10 +51,6 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& ! nodata - number of BUFR IASI-NG observations retained for further processing ! nobs - array of observations on each subdomain for each processor ! -! attributes: -! language: f90 -! machine: ibm RS/6000 SP -! !$$$ ! Use modules use kinds, only: r_kind,r_double,i_kind @@ -157,7 +154,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& logical :: quiet,cloud_info integer(i_kind) :: ifov, ifor, istep, ipos, instr, iscn, ioff, sensorindex_iasing - integer(i_kind) :: i, j, l, iskip, ifovn, bad_line, ksatid, kidsat, llll + integer(i_kind) :: i, j, l, iskip, ifovn, ksatid, kidsat, llll integer(i_kind) :: nreal, isflg integer(i_kind) :: itx, k, nele, itt, n integer(i_kind) :: iexponent,maxinfo, bufr_nchan, dval_info @@ -213,8 +210,6 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& ndata = 0 nodata = 0 - bad_line=-1 - if (nst_gsi > 0 ) then call gsi_nstcoupler_skindepth(obstype, zob) ! get penetration depth (zob) for the obstype endif @@ -447,13 +442,6 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& if ( linele(3) /= zero) cycle read_loop ! problem with profile (INGGQF) -! Zenith angle/scan spot mismatch, reject entire line - if ( bad_line == nint(linele(2))) then - cycle read_loop - else - bad_line = -1 - endif - ifov = nint(linele(1)) ! field of view if ( ifov < 1 .or. ifov > 224 ) then ! field of view out of range write(6,*)'READ_IASI-NG: ### ERROR IN READING ', senname, ' BUFR DATA:', & From 745f01d94048ea13fa1e20ef9d6b86b580173ee0 Mon Sep 17 00:00:00 2001 From: wx20jjung Date: Tue, 5 Nov 2024 17:42:50 +0000 Subject: [PATCH 13/16] Added comments around some of the bufr read statements. --- src/gsi/read_iasing.f90 | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/gsi/read_iasing.f90 b/src/gsi/read_iasing.f90 index c20de10524..3443e76557 100644 --- a/src/gsi/read_iasing.f90 +++ b/src/gsi/read_iasing.f90 @@ -419,11 +419,12 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& read_loop: do while (ireadsb(lnbufr)==0) ! Get the size of the channels and radiance (allchan) array +! This is a delayed replication. crchn_reps is the number of IASI-NG replications (channel and radiance) call ufbint(lnbufr,crchn_reps,1,1,iret,'(I1CRSQ)') bufr_nchan = int(crchn_reps) bufr_size = size(temperature,1) - if ( bufr_size /= bufr_nchan ) then ! Re-allocation if number of channels has changed + if ( bufr_size /= bufr_nchan ) then ! Re-allocation if number of channels has changed in the delayed replication ! Allocate the arrays needed for the channel and radiance array deallocate(temperature,allchan,bufr_chan_test,scalef) allocate(temperature(bufr_nchan)) ! dependent on # of channels in the bufr file @@ -638,10 +639,12 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& endif end do -! Read IASI-NG channel number(CHNM) and radiance (SCRA) +! Read IASI-NG channel number(CHNM) and radiance (SCRA). call ufbseq(lnbufr,allchan,2,bufr_nchan,iret,'I1CRSQ') jstart=1 scalef=one + +! Determine the scaling factor (scalef(i)) for the radiance (allchan(2,i)) of each channel (allchan(1,i)) do i=1,bufr_nchan scaleloop: do j=jstart,4 if(allchan(1,i) >= cscale(1,j) .and. allchan(1,i) <= cscale(2,j))then @@ -729,6 +732,7 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& if(.not. iuse)cycle read_loop ! Read the imager cluster information for the Cloud and Aerosol Detection Software. +! This is a nested fixed replication read of the metimage data. ! Only channels 18 and 19 are used. if ( iasing_cads ) then From c4afc9987e2fb165ffbedfe2c2f7834cc8fa5a63 Mon Sep 17 00:00:00 2001 From: wx20jjung Date: Tue, 5 Nov 2024 22:07:02 +0000 Subject: [PATCH 14/16] Added comments --- src/gsi/read_iasing.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/gsi/read_iasing.f90 b/src/gsi/read_iasing.f90 index 3443e76557..7de2f67bb9 100644 --- a/src/gsi/read_iasing.f90 +++ b/src/gsi/read_iasing.f90 @@ -298,7 +298,8 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& ! find IASI-NG sensorindex sensorindex_iasing = 0 if ( sc(1)%sensor_id(1:7) == 'iasi-ng' .or. sc(1)%sensor_id == '999') then -!JAJ if ( sc(1)%sensor_id(1:7) == 'iasi-ng' ) then +!TODO This is the correct code once the CRTM spectral coefficient is fixed. +!TODO JAJ if ( sc(1)%sensor_id(1:7) == 'iasi-ng' ) then sensorindex_iasing = 1 else write(6,*)'READ_IASI-NG: ***ERROR*** sensorindex_iasi-ng not set NO IASI-NG DATA USED' @@ -560,7 +561,6 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& ! | 9 10 11 12 | ! | 13 14 15 16 | -! Add this code when field-of-view spans 1-224 JAJ. ! There are 14 fields-of-regard in each scan line. So, there are 56 unique positions in the scan line. ! To determine the scan position: ifor = (ifov-1) / 16 ! Determine field-of-regard From 5db524cc941157d218b4725efbbceadd0ced7b74 Mon Sep 17 00:00:00 2001 From: "jim.jung" Date: Wed, 6 Nov 2024 20:01:58 +0000 Subject: [PATCH 15/16] Added/changed comments --- src/gsi/crtm_interface.f90 | 8 ++++++-- src/gsi/read_iasing.f90 | 7 +++++-- 2 files changed, 11 insertions(+), 4 deletions(-) diff --git a/src/gsi/crtm_interface.f90 b/src/gsi/crtm_interface.f90 index 2f2764b827..db8ae6298b 100644 --- a/src/gsi/crtm_interface.f90 +++ b/src/gsi/crtm_interface.f90 @@ -667,8 +667,12 @@ subroutine init_crtm(init_pass,mype_diaghdr,mype,nchanl,nreal,isis,obstype,radmo error_status = crtm_channelinfo_subset(channelinfo(1), & channel_subset = nuchan(subset_start:subset_end)) -!TODO JAJ the correct one else if (channelinfo(1)%sensor_id(1:7) == 'iasi-ng' .AND. isis(1:7) == 'iasi-ng') then -else if (channelinfo(1)%sensor_id(1:3) == '999' .AND. isis(1:7) == 'iasi-ng') then ! TODO To be removed. +! TODO The CRTM spectral coefficient files have the instrument name in the beginning of the file. The current iasi-ng coefficient +! TODO file contains '999' instead of the instrument name. When the final coefficient file is built, it will have 'iasi-ng'. +! TODO else if (channelinfo(1)%sensor_id(1:7) == 'iasi-ng' .AND. isis(1:7) == 'iasi-ng') then +! TODO when this file exists, use the above line. +else if (channelinfo(1)%sensor_id(1:3) == '999' .AND. isis(1:7) == 'iasi-ng') then +! TODO and remove the above line. sensorindex = 1 subset_start = 0 subset_end = 0 diff --git a/src/gsi/read_iasing.f90 b/src/gsi/read_iasing.f90 index 7de2f67bb9..4183c59843 100644 --- a/src/gsi/read_iasing.f90 +++ b/src/gsi/read_iasing.f90 @@ -297,9 +297,12 @@ subroutine read_iasing(mype,val_iasing,ithin,isfcalc,rmesh,jsatid,gstime,& ! find IASI-NG sensorindex sensorindex_iasing = 0 +!TODO The final iasi-ng spectral coefficient file is not available yet. The current file has '999' in place of iasi-ng. +!TODO The line below is a temporary fix to work around this problem in the CRTM spectral coefficient file. if ( sc(1)%sensor_id(1:7) == 'iasi-ng' .or. sc(1)%sensor_id == '999') then -!TODO This is the correct code once the CRTM spectral coefficient is fixed. -!TODO JAJ if ( sc(1)%sensor_id(1:7) == 'iasi-ng' ) then +!TODO Remove the above line when the correct CRTM file is available +!TODO if ( sc(1)%sensor_id(1:7) == 'iasi-ng' ) then +!TODO Use the above line once the CRTM spectral coefficient file is fixed sensorindex_iasing = 1 else write(6,*)'READ_IASI-NG: ***ERROR*** sensorindex_iasi-ng not set NO IASI-NG DATA USED' From 493e8e65954e48db7c183e25d3df60b42a3cf68a Mon Sep 17 00:00:00 2001 From: James Date: Thu, 7 Nov 2024 18:15:05 +0000 Subject: [PATCH 16/16] Changes to Jet node/task configuration to run ctests. --- regression/regression_param.sh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/regression/regression_param.sh b/regression/regression_param.sh index 34d5964d87..c59a4285c3 100755 --- a/regression/regression_param.sh +++ b/regression/regression_param.sh @@ -101,8 +101,8 @@ case $regtest in topts[1]="0:05:00" ; popts[1]="40/3/" ; ropts[1]="/1" topts[2]="0:05:00" ; popts[2]="40/5/" ; ropts[2]="/2" elif [[ "$machine" = "Jet" ]]; then - topts[1]="0:15:00" ; popts[1]="5/4/" ; ropts[1]="/1" - topts[2]="0:15:00" ; popts[2]="10/4/" ; ropts[2]="/1" + topts[1]="0:15:00" ; popts[1]="40/3/" ; ropts[1]="/1" + topts[2]="0:15:00" ; popts[2]="40/5/" ; ropts[2]="/1" elif [[ "$machine" = "Gaea" ]]; then topts[1]="0:15:00" ; popts[1]="64/1/" ; ropts[1]="/1" topts[2]="0:15:00" ; popts[2]="128/2/" ; ropts[2]="/1"