From 2b1d706fa4af91ba1e9229fef217e2fba6e5867d Mon Sep 17 00:00:00 2001 From: James Jung Date: Fri, 20 Sep 2024 11:30:36 -0400 Subject: [PATCH] Iasi debug fix (#790) --- src/enkf/enkf_obs_sensitivity.f90 | 2 +- src/enkf/innovstats.f90 | 2 +- src/gsi/combine_radobs.f90 | 10 +++++----- src/gsi/radinfo.f90 | 10 +++++----- src/gsi/read_avhrr_navy.f90 | 1 + src/gsi/read_bufrtovs.f90 | 2 +- src/gsi/read_cris.f90 | 2 +- src/gsi/read_gmi.f90 | 1 - src/gsi/read_iasi.f90 | 21 +++++++++------------ src/gsi/read_saphir.f90 | 3 +-- src/gsi/statsrad.f90 | 2 +- 11 files changed, 26 insertions(+), 30 deletions(-) diff --git a/src/enkf/enkf_obs_sensitivity.f90 b/src/enkf/enkf_obs_sensitivity.f90 index 72296d5934..1b95b45a7a 100644 --- a/src/enkf/enkf_obs_sensitivity.f90 +++ b/src/enkf/enkf_obs_sensitivity.f90 @@ -547,7 +547,7 @@ subroutine print_ob_sens if(nob_sat(nchan) > 0) then rate_sat(nchan,1:3) = rate_sat(nchan,1:3) & & / real(nob_sat(nchan),r_kind) * 100._r_kind - write(*,'(a20,i5,i7,3(1x,e12.5),3(1x,f7.2))') & + write(*,'(a20,1x,i5,i7,3(1x,e12.5),3(1x,f7.2))') & & trim(adjustl(nusis(nchan))), & & nuchan(nchan),nob_sat(nchan),sumsense_sat(nchan,1:3), & & rate_sat(nchan,1:3) diff --git a/src/enkf/innovstats.f90 b/src/enkf/innovstats.f90 index 853532c9b9..0abd12a919 100644 --- a/src/enkf/innovstats.f90 +++ b/src/enkf/innovstats.f90 @@ -278,7 +278,7 @@ subroutine print_innovstats(obfit,obsprd) sqrt(sumerr_sat(nchan)) end if end do -9805 format(a20,i4,1x,i5,5(1x,e10.3)) +9805 format(a20,1x,i5,1x,i5,5(1x,e10.3)) end if !nobs_sat>0 end subroutine print_innovstats diff --git a/src/gsi/combine_radobs.f90 b/src/gsi/combine_radobs.f90 index 7692bdef3b..b680c35648 100644 --- a/src/gsi/combine_radobs.f90 +++ b/src/gsi/combine_radobs.f90 @@ -25,11 +25,10 @@ subroutine combine_radobs(mype_sub,mype_root,& ! data_all - observation data array ! data_crit- array containing observation "best scores" ! nread - task specific number of obesrvations read from data file -! ndata - task specific number of observations keep for assimilation ! ! output argument list: ! nread - total number of observations read from data file (mype_root) -! ndata - total number of observations keep for assimilation (mype_root) +! ndata - total number of observation profiles kept for assimilation in the thinning box (mype_root) ! data_all - merged observation data array (mype_root) ! data_crit- merged array containing observation "best scores" (mype_root) ! @@ -50,7 +49,8 @@ subroutine combine_radobs(mype_sub,mype_root,& integer(i_kind) ,intent(in ) :: npe_sub,itxmax integer(i_kind) ,intent(in ) :: nele integer(i_kind) ,intent(in ) :: mpi_comm_sub - integer(i_kind) ,intent(inout) :: nread,ndata + integer(i_kind) ,intent(inout) :: nread + integer(i_kind) ,intent( out) :: ndata integer(i_kind),dimension(itxmax) ,intent(in ) :: nrec real(r_kind),dimension(itxmax) ,intent(inout) :: data_crit real(r_kind),dimension(nele,itxmax),intent(inout) :: data_all @@ -74,7 +74,7 @@ subroutine combine_radobs(mype_sub,mype_root,& nread=0 if (mype_sub==mype_root) nread = ncounts1 - if (ncounts1 == 0)return + if (ncounts1 <= 0)return ! Allocate arrays to hold data @@ -83,7 +83,7 @@ subroutine combine_radobs(mype_sub,mype_root,& ! is only needed on task mype_root call mpi_allreduce(data_crit,data_crit_min,itxmax,mpi_rtype,mpi_min,mpi_comm_sub,ierror) - allocate(nloc(min(ncounts1,itxmax)),icrit(min(ncounts1,itxmax))) + allocate(nloc(itxmax),icrit(itxmax)) icrit=1e9 ndata=0 ndata1=0 diff --git a/src/gsi/radinfo.f90 b/src/gsi/radinfo.f90 index 8bb496c645..e7624fab70 100644 --- a/src/gsi/radinfo.f90 +++ b/src/gsi/radinfo.f90 @@ -805,7 +805,7 @@ subroutine radinfo_read end do close(lunin) 100 format(a1,a120) -110 format(i4,1x,a20,' chan= ',i5, & +110 format(i5,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= ',i5, & @@ -1135,7 +1135,7 @@ subroutine radinfo_read nusis(j),nuchan(j),' not found in satbias_in file - set to zero ' endif end do -140 format(i4,1x,a20,12f12.6) +140 format(i5,1x,a20,12f12.6) endif @@ -1687,7 +1687,6 @@ subroutine init_predx integer(i_kind),parameter:: lntemp = 51 integer(i_kind),parameter:: nthreshold = 100 - integer(i_kind),parameter:: maxchn = 3000 integer(i_kind),parameter:: maxdat = 100 real(r_kind), parameter:: atiny = 1.0e-10_r_kind @@ -1712,7 +1711,7 @@ subroutine init_predx integer(i_kind):: np,new_chan,nc integer(i_kind):: counttmp, jjstart, sensor_start, sensor_end integer(i_kind):: radedge_min, radedge_max - integer(i_kind),dimension(maxchn):: ich + integer(i_kind),allocatable,dimension(:):: ich integer(i_kind),dimension(maxdat):: ipoint real(r_kind):: bias,scan,errinv,rnad @@ -1814,6 +1813,7 @@ subroutine init_predx mype, trim(fdiag_rad), header_fix%idate satsens = header_fix%isis n_chan = header_fix%nchan + allocate(ich(n_chan)) ! Check for consistency between specified and retrieved satellite id ! after first sorting out some historical naming conventions @@ -2063,7 +2063,7 @@ subroutine init_predx if ( nuchan(jj) == header_chan(j)%nuchan ) then jjstart = jj + 1 write(lntemp,220) jj,tlapmean(jj),tsum_tlapmean(jj),count_tlapmean(jj) -220 format(I5,1x,2e15.6,1x,I5) +220 format(I5,1x,2e15.6,1x,I6) cycle loop_c endif end do diff --git a/src/gsi/read_avhrr_navy.f90 b/src/gsi/read_avhrr_navy.f90 index dd5a64083a..b21e653730 100644 --- a/src/gsi/read_avhrr_navy.f90 +++ b/src/gsi/read_avhrr_navy.f90 @@ -255,6 +255,7 @@ subroutine read_avhrr_navy(mype,val_avhrr,ithin,rmesh,jsatid,& next=0 ! Read BUFR Navy data + nrec = 999999 irec=0 read_msg: do while (ireadmg(lnbufr,subset,idate) >= 0) irec=irec+1 diff --git a/src/gsi/read_bufrtovs.f90 b/src/gsi/read_bufrtovs.f90 index 2fc14b5cdf..9bca8f5a7c 100644 --- a/src/gsi/read_bufrtovs.f90 +++ b/src/gsi/read_bufrtovs.f90 @@ -490,7 +490,7 @@ subroutine read_bufrtovs(mype,val_tovs,ithin,isfcalc,& hdr2b ='SAZA SOZA BEARAZ SOLAZI' allocate(data_all(nele,itxmax),data1b8(nchanl),data1b4(nchanl),nrec(itxmax)) - + nrec = 999999 next=0 irec=0 ! Big loop over standard data feed and possible ears/db data diff --git a/src/gsi/read_cris.f90 b/src/gsi/read_cris.f90 index 84288a7f04..2a899cecd6 100644 --- a/src/gsi/read_cris.f90 +++ b/src/gsi/read_cris.f90 @@ -455,7 +455,7 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& ! Big loop to read data file next=0 irec=0 - nrec = 99999 + nrec = 999999 ! Big loop over standard data feed and possible rars/db data ! llll=1 is normal feed, llll=2 RARS data, llll=3 DB/UW data) ears_db_loop: do llll= 1, 3 diff --git a/src/gsi/read_gmi.f90 b/src/gsi/read_gmi.f90 index 6ad4d829a3..d02520d7e0 100644 --- a/src/gsi/read_gmi.f90 +++ b/src/gsi/read_gmi.f90 @@ -346,7 +346,6 @@ subroutine read_gmi(mype,val_gmi,ithin,rmesh,jsatid,gstime,& next=0 irec=0 iobs=1 - nrec=999999 read_subset: do while(ireadmg(lnbufr,subset,idate)>=0) ! GMI scans irec=irec+1 diff --git a/src/gsi/read_iasi.f90 b/src/gsi/read_iasi.f90 index 75dc50bb76..94e1577798 100644 --- a/src/gsi/read_iasi.f90 +++ b/src/gsi/read_iasi.f90 @@ -759,17 +759,16 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& !$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_iasi,sc_chan,radiance,temperature(bufr_chan)) - else - temperature(bufr_chan) = tbmin - endif - end if + 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_iasi,sc_chan,radiance,temperature(bufr_chan)) + else + temperature(bufr_chan) = tbmin + endif end do channel_loop ! Check for reasonable temperature values @@ -950,10 +949,8 @@ subroutine read_iasi(mype,val_iasi,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 diff --git a/src/gsi/read_saphir.f90 b/src/gsi/read_saphir.f90 index 06e992b03d..dd421fd679 100644 --- a/src/gsi/read_saphir.f90 +++ b/src/gsi/read_saphir.f90 @@ -110,7 +110,7 @@ subroutine read_saphir(mype,val_tovs,ithin,isfcalc,& character(8) :: subset character(80) :: hdr1b,hdr2b - integer(i_kind) :: ireadsb,ireadmg,irec + integer(i_kind) :: ireadsb,ireadmg integer(i_kind) :: i,j,k,ntest,iob integer(i_kind) :: iret,idate,nchanl,n,idomsfc(1) integer(i_kind) :: kidsat,maxinfo @@ -293,7 +293,6 @@ subroutine read_saphir(mype,val_tovs,ithin,isfcalc,& ! hdr2b ='AGIND SOZA BEARAZ SOLAZI' ! AGIND instead of SAZA ! Loop to read bufr file - irec=0 read_subset: do while(ireadmg(lnbufr,subset,idate)>=0 .AND. iob < maxobs) read_loop: do while (ireadsb(lnbufr)==0 .and. iob < maxobs) diff --git a/src/gsi/statsrad.f90 b/src/gsi/statsrad.f90 index dfedcc92ab..986ab4e226 100644 --- a/src/gsi/statsrad.f90 +++ b/src/gsi/statsrad.f90 @@ -161,7 +161,7 @@ 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,i6,1x,a20,2i7,1x,f10.3,1x,6(f11.7,1x)) +1102 format(1x,i6,i6,1x,a20,2i7,1x,f10.3,1x,6(f11.7,1x)) 1109 format(t5,'it',t13,'satellite',t23,'instrument',t40, & '# read',t53,'# keep',t65,'# assim',& t75,'penalty',t88,'qcpnlty',t104,'cpen',t115,'qccpen')