Iasi debug fix

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)

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
 

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

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(i6,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

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

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

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

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

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)
 

src/gsi/read_seviri.f90

@@ -103,7 +103,7 @@ subroutine read_seviri(mype,val_sev,ithin,rmesh,jsatid,&
   integer(i_kind) nchanl,ilath,ilonh,ilzah,iszah,irec,next
   integer(i_kind) nmind,lnbufr,idate,ilat,ilon,nhdr,nchn,ncld,nbrst,jj
   integer(i_kind) ireadmg,ireadsb,iret,nreal,nele,itt
-  integer(i_kind) itx,i,k,isflg,kidsat,n,iscan,idomsfc
+  integer(i_kind) itx,i,k,isflg,kidsat,n,iscan,idomsfc,number_profiles
   integer(i_kind) idate5(5),maxinfo
   integer(i_kind),allocatable,dimension(:)::nrec
 
@@ -528,8 +528,10 @@ subroutine read_seviri(mype,val_sev,ithin,rmesh,jsatid,&
   call closbf(lnbufr)
   close(lnbufr)
 
+  number_profiles = count(nrec(:) /= 999999,dim=1)
+
   call combine_radobs(mype_sub,mype_root,npe_sub,mpi_comm_sub,&
-     nele,itxmax,nread,ndata,data_all,score_crit,nrec)
+     nele,itxmax,nread,number_profiles,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.

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')