Changes to code based on shoyokota changes #3

src/gsi/convthin.f90

@@ -210,9 +210,9 @@ subroutine createfore
 !
 !$$$
     integer i,j
-    allocate(icount(itxmax,nlevp))
-    allocate(ibest_obs(itxmax,nlevp))
-    allocate(score_crit(itxmax,nlevp))
+    allocate(icount_fore(itxmax,nlevp))
+    allocate(ibest_obs_fore(itxmax,nlevp))
+    allocate(score_crit_fore(itxmax,nlevp))
 
     do j=1,nlevp
        do i=1,itxmax
@@ -227,7 +227,7 @@ end subroutine createfore
   subroutine createaft
 !$$$  subprogram documentation block
 !                .      .    .                                       .
-! subprogram:    createnormal
+! subprogram:    createaft
 !     prgmmr:    derber     org: np23                date: 2023-10-20
 !
 ! abstract:  This routine creates and initializes arrays for aft thinning

src/gsi/convthin_time.f90

@@ -189,7 +189,7 @@ end subroutine createnormal_tm
   subroutine createfore_tm
 !$$$  subprogram documentation block
 !                .      .    .                                       .
-! subprogram:    createnormal
+! subprogram:    createfore_tm
 !     prgmmr:    derber     org: np23                date: 2023-10-20
 !
 ! abstract:  This routine creates and initializes arrays for fore thinning

src/gsi/ensctl2state.f90

@@ -49,7 +49,7 @@ module ensctl2state_mod
 logical :: ls_w,ls_dw
 
 logical :: lc_sf,lc_vp,lc_ps,lc_t,lc_rh,lc_cw
-logical :: lc_w,lc_dw,wdw_exist
+logical :: lc_w,lc_dw
 
 logical :: do_getuv,do_tv_to_tsen,do_normal_rh_to_q,do_getprs,lstrong_bk_vars
 logical :: do_q_copy
@@ -233,20 +233,15 @@ subroutine ensctl2state(xhat,mval,eval)
 
 !  Get pointers to required state variables and copy
    call gsi_bundlegetpointer (eval(jj),'sst' ,sv_sst, istatus)
-   if(ls_w)then
+   call gsi_bundlegetvar ( wbundle_c, 'sst', sv_sst, istatus )
+   if(ls_w .and. lc_w)then
      call gsi_bundlegetpointer (eval(jj),'w' ,sv_w, istatus)
-     if(ls_dw)then
+     call gsi_bundlegetvar ( wbundle_c, 'w' , sv_w,  istatus )
+     if(ls_dw .and. lc_dw)then
         call gsi_bundlegetpointer (eval(jj),'dw' ,sv_dw, istatus)
+        call gsi_bundlegetvar ( wbundle_c, 'dw' , sv_dw,  istatus )
      end if
    end if
-!  Copy variables
-   call gsi_bundlegetvar ( wbundle_c, 'sst', sv_sst, istatus )
-   if(lc_w)then
-      call gsi_bundlegetvar ( wbundle_c, 'w' , sv_w,  istatus )
-      if(lc_dw)then
-         call gsi_bundlegetvar ( wbundle_c, 'dw' , sv_dw,  istatus )
-      end if
-   end if
 
 !  Get the ozone vector if it is defined
    if(idozone > 0) then
@@ -355,7 +350,6 @@ subroutine ensctl2state_set(xhat,eval)
 do_cw_to_hydro = lc_cw .and. ls_ql .and. ls_qi
 do_cw_to_hydro_hwrf = lc_cw.and.ls_ql.and.ls_qi.and.ls_qr.and.ls_qs.and.ls_qg.and.ls_qh
 
-wdw_exist = lc_w.and.lc_dw.and.ls_w.and.ls_dw
 
 idozone=getindex(cvars3d,"oz")
 
@@ -496,10 +490,10 @@ subroutine ensctl2state_ad(eval,mval,grad)
 
    call gsi_bundlegetpointer (eval(jj),'sst' ,rv_sst, istatus)
    call gsi_bundleputvar ( wbundle_c, 'sst', rv_sst, istatus )
-   if(wdw_exist)then
+   if(lc_w .and. ls_w)then
      call gsi_bundlegetpointer (eval(jj),'w' ,rv_w, istatus)
      call gsi_bundleputvar ( wbundle_c, 'w', rv_w, istatus )
-     if(ls_dw)then
+     if(ls_dw .and. lc_dw)then
        call gsi_bundlegetpointer (eval(jj),'dw' ,rv_dw, istatus)
        call gsi_bundleputvar ( wbundle_c, 'dw', rv_dw, istatus )
      end if

src/gsi/read_dbz_nc.f90

@@ -150,7 +150,7 @@ subroutine read_dbz_nc(nread,ndata,nodata,infile,lunout,obstype,sis,hgtl_full,no
                   dlat,dlon,thiserr,thislon,thislat, &
                   timeb
   real(r_kind) :: radartwindow
-  real(r_kind) :: rmins_an
+  real(r_kind) :: rmins_an,usage
   real(r_kind),allocatable,dimension(:,:):: cdata_all
   real(r_double) rstation_id
   logical, allocatable,dimension(:)     :: rusage,rthin
@@ -228,8 +228,6 @@ subroutine read_dbz_nc(nread,ndata,nodata,infile,lunout,obstype,sis,hgtl_full,no
   icntpnt=0
   zflag=0
 
-  rusage = .true.
-  rthin = .false.
   use_all=.true.
   if (ithin > 0) then
      write(6,*)'READ_RADAR_DBZ: ithin,rmesh :',ithin,rmesh
@@ -431,7 +429,8 @@ subroutine read_dbz_nc(nread,ndata,nodata,infile,lunout,obstype,sis,hgtl_full,no
 
         nread = nread + 1
 
-        if(icuse(ikx) < zero)rusage(ndata+1)=.false.
+        usage=zero
+        if(icuse(ikx) < zero)usage=r100
  !####################       Data thinning       ###################
         icntpnt=icntpnt+1
         if(icntpnt>maxobs) exit
@@ -517,6 +516,7 @@ subroutine read_dbz_nc(nread,ndata,nodata,infile,lunout,obstype,sis,hgtl_full,no
          cdata_all(17,iout)= dbznoise                      ! noise threshold for reflectivity (dBZ)
          cdata_all(18,iout)= imissing2nopcp                !=0, normal 
                                                            !=1,  !values !converted !from !missing !values 
+         if(usage >= r100)rusage(ndata)=.false.
 
          if(doradaroneob .and. (cdata_all(5,iout) > -99.0_r_kind) ) exit ILOOP
 
@@ -592,10 +592,10 @@ subroutine read_dbz_nc(nread,ndata,nodata,infile,lunout,obstype,sis,hgtl_full,no
 
   !---------------DEALLOCATE ARRAYS-------------!
 
-  deallocate(cdata_all,rusage,rthin)
 else  !fileopen
   write(6,*) 'READ_dBZ: ERROR OPENING RADAR REFLECTIVITY FILE: ',trim(infile),' IOSTAT ERROR: ',ierror, ' SKIPPING...'
 end if fileopen
+deallocate(cdata_all,rusage,rthin)
 
 
 end subroutine read_dbz_nc

src/gsi/read_dbz_netcdf.f90

@@ -605,20 +605,20 @@ subroutine read_dbz_mrms_netcdf(nread,ndata,nodata,infile,obstype,lunout,sis,nob
 
   !---------------DEALLOCATE ARRAYS-------------!
 
-  deallocate(cdata_all)
-  do v=1,nvol
-     do k=1,nelv
-        deallocate(strct_in_dbz(v,k)%azim)
-        deallocate(strct_in_dbz(v,k)%field)
-     end do
-  end do
-  deallocate(strct_in_dbz)
-  deallocate(obdata_nc)
-  deallocate(beamwidth_nc,azimspacing_nc,gatewidth_nc)
 
  else  !fileopen
   write(6,*) 'READ_dBZ: ERROR OPENING RADAR REFLECTIVITY FILE: ',trim(infile),' IOSTAT ERROR: ',ierror, ' SKIPPING...'
  end if fileopen
+ deallocate(cdata_all)
+ do v=1,nvol
+    do k=1,nelv
+       deallocate(strct_in_dbz(v,k)%azim)
+       deallocate(strct_in_dbz(v,k)%field)
+    end do
+ end do
+ deallocate(strct_in_dbz)
+ deallocate(obdata_nc,azimuth_nc)
+ deallocate(beamwidth_nc,azimspacing_nc,gatewidth_nc)
 
 end subroutine read_dbz_mrms_netcdf
 

src/gsi/read_fl_hdob.f90

@@ -672,10 +672,7 @@ subroutine read_fl_hdob(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,si
 !          temperature will be used
            if (ltob) then
               qcm = t_qm 
-              if (qcm > 3) then
-                 usage = r100  
-                 if(pmot >= 2 .and. usage >= r100)rusage(ndata+1)=.false.
-              end if
+              if (qcm > 3) usage = r100  
               call ufbint(lunin,obstmp,2,1,nlv,tmpstr)
               call ufbint(lunin,obsmst,3,1,nlv,mststr)
               tob  = obstmp(2,1)  ! airs temperature [K] 
@@ -755,10 +752,7 @@ subroutine read_fl_hdob(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,si
 !          related QC mark and then calculate specific humidity
            if (lqob) then
               qcm = q_qm
-              if (qcm > 3) then
-                 usage = r100 
-                 if(pmot >= 2 .and. usage >= r100)rusage(ndata+1)=.false.
-              end if
+              if (qcm > 3) usage = r100 
               call ufbint(lunin,obstmp,2,1,nlv,tmpstr)
               call ufbint(lunin,obsmst,3,1,nlv,mststr)
               tob  = obstmp(2,1)  ! dry airs temperature [K] 
@@ -831,10 +825,7 @@ subroutine read_fl_hdob(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,si
 !          Convert wind direction and spped to u and v wind components
            if (luvob) then
               qcm = uv_qm
-              if (qcm > 3) then
-                 usage=r100 
-                 if(pmot >= 2 .and. usage >= r100)rusage(ndata+1)=.false.
-              end if
+              if (qcm > 3) usage=r100 
               call ufbint(lunin,obswnd,4,1,nlv,wndstr)
               wdir = obswnd(2,1)     ! degree true  
               wspd = obswnd(3,1)     ! m/s 
@@ -846,10 +837,7 @@ subroutine read_fl_hdob(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,si
 !          Read surface wind speed [m/s] and total rain rate [mm/hr] from SFMR 
            if (lspdob) then
               qcm = wspd_qm
-              if (qcm > 3) then
-                 usage=r100 
-                 if(pmot >= 2 .and. usage >= r100)rusage(ndata+1)=.false.
-              end if
+              if (qcm > 3) usage=r100 
               call ufbint(lunin,obsfmr,2,1,nlv,sfmrstr)
 !             print*, 'PKSWSP =  ', obstype,obsfmr(1,1)
 !             print*, 'TRRP   =  ', obstype,obsfmr(2,1)
@@ -1164,6 +1152,7 @@ subroutine read_fl_hdob(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,si
               cdata_all(22,iout)=r_prvstg(1,1)          !  provider name 
               cdata_all(23,iout)=r_sprvstg(1,1)         !  subprovider name 
            endif 
+           if(usage >= r100)rusage(ndata)=.false.
 
         end do loop_readsb2
      end do loop_msg2
@@ -1194,7 +1183,7 @@ subroutine read_fl_hdob(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,si
 !       end do
 !       write(6,*) ' fl ',trim(ioctype(nc)),ictype(nc),icsubtype(nc),numall,&
 !              numrem,numqc,numthin
-!   If thinned data set quality mark to 16
+!   If thinned data set quality mark to 14
         if (ithin > 0 .and. ithin <5) then
           do i=1,nxdata
              if(rthin(i))cdata_all(iqm,i)=14
@@ -1232,7 +1221,7 @@ subroutine read_fl_hdob(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,si
 !    Write header record and data to output file for further processing
 !     deallocate(etabl)
 
-     call count_obs(ndata,nreal,ilat,ilon,cdata_all,nobs)
+     call count_obs(ndata,nreal,ilat,ilon,cdata_all(1,1:ndata),nobs)
      write(lunout) obstype,sis,nreal,nchanl,ilat,ilon
      write(lunout) ((cdata_all(k,i),k=1,nreal),i=1,ndata)
      deallocate(cdata_all,rusage,rthin)

src/gsi/read_gfs_ozone_for_regional.f90

@@ -413,7 +413,7 @@ subroutine read_gfs_ozone_for_regional
         write(6,*)'READ_GFS_OZONE_FOR_REGIONAL:  ***ERROR*** INVALID value for nvcoord=',sighead%nvcoord,filename
         call stop2(85)
      endif
-   else if ( use_gfs_ncio ) then 
+  else if ( use_gfs_ncio ) then 
      if (gfshead%nvcoord == 1) then
         do k=1,nsig_gfs+1
            bk5(k) = gfsheadv%vcoord(k,1)
@@ -455,9 +455,8 @@ subroutine read_gfs_ozone_for_regional
         write(6,*)'GET_GEFS_FOR_REGIONAL:  ***ERROR*** INVALID value for nvcoord=',nvcoord
         call stop2(85)
      endif
-     deallocate(vcoord)
+     deallocate(vcoord,nems_vcoord)
   end if
-  deallocate(nems_vcoord)
 
 ! Load reference temperature array (used by general coordinate)
   do k=1,nsig_gfs
@@ -507,8 +506,6 @@ subroutine read_gfs_ozone_for_regional
 !  also want to set up regional grid structure variable grd_mix, which still has number of
 !   vertical levels set to nsig_gfs, but horizontal dimensions set to regional domain.
 
-  num_fields=2*nsig_gfs
-  vector=.false.
   call general_sub2grid_create_info(grd_mix,inner_vars,nlat,nlon,nsig_gfs,num_fields,regional,vector)
   deallocate(vector)
 

src/gsi/read_goesimgr_skycover.f90

@@ -340,7 +340,6 @@ subroutine  read_goesimgr_skycover(nread,ndata,nodata,infile,obstype,lunout,gsti
       !-  Set usage variable
       usage = 0 
       if(iuse <= 0)usage=r100
-      if(usage >= r100 .and. pmot >= 2)rusage(ndata+1)=.false.
 
       ! Get information from surface file necessary for conventional data here
       call deter_sfc2(dlat_earth,dlon_earth,t4dv,idomsfc,tsavg,ff10,sfcr,zz)
@@ -378,6 +377,7 @@ subroutine  read_goesimgr_skycover(nread,ndata,nodata,infile,obstype,lunout,gsti
       cdata_all(18,iout)=zz                     !  terrain height at ob location
       cdata_all(19,iout)=r_prvstg(1,1)          !  provider name
       cdata_all(20,iout)=r_sprvstg(1,1)         !  subprovider name
+      if(usage >=r100)rusage(ndata)=.false.
 
      enddo  loop_readsb
 
@@ -441,7 +441,7 @@ subroutine  read_goesimgr_skycover(nread,ndata,nodata,infile,obstype,lunout,gsti
 
 ! Write header record and data to output file for further processing
 
-  call count_obs(ndata,nreal,ilat,ilon,cdata_all,nobs)
+  call count_obs(ndata,nreal,ilat,ilon,cdata_all(1,1:ndata),nobs)
   write(lunout) obstype,sis,nreal,nchanl,ilat,ilon,ndata
   write(lunout) ((cdata_all(k,i),k=1,nreal),i=1,ndata)
 

src/gsi/read_prepbufr.f90

@@ -2889,9 +2889,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
                                   hig_cldamt,hig_cldamt_qc,tcamt,lcbas,tcamt_qc,lcbas_qc,ceiling,stnelev)
 
                     if(lcbas_qc >= 8) usage=100._r_kind
-                    if(pmot >=2 .and. usage >= 100.0_r_kind)then
-                       rusage(iout)=.false.
-                    end if
+                    if(usage >= 100.0_r_kind)rusage(iout)=.false.
                     lcbas_oe=4500.0_r_kind
                     if(lcbas_qc==3) lcbas_oe=lcbas_oe*1.25_r_kind
                     if(lcbas_qc==4) lcbas_oe=lcbas_oe*1.5_r_kind

src/gsi/read_radar.f90

@@ -862,12 +862,11 @@ subroutine read_radar(nread,ndata,nodata,infile,lunout,obstype,twind,sis,hgtl_fu
            if(ncnumgrp(ikx) > 0 )then                     ! cross validation on
               if(mod(ndata,ncnumgrp(ikx))== ncgroup(ikx)-1)usage=ncmiter(ikx)
            end if
-           if(pmot >=2 .and. usage >= 100._r_kind) rusage(ndata)=.false.
+           if(usage >= 100._r_kind) rusage(ndata)=.false.
 
            call deter_sfc2(dlat_earth,dlon_earth,t4dv,idomsfc,skint,ff10,sfcr)
 
            LEVEL_TWO_READ: if(loop==0 .and. sis=='l2rw') then      
-              write(6,*) ' radar1 ',height
               cdata(1) = error             ! wind obs error (m/s)
               cdata(2) = dlon              ! grid relative longitude
               cdata(3) = dlat              ! grid relative latitude
@@ -1302,7 +1301,6 @@ subroutine read_radar(nread,ndata,nodata,infile,lunout,obstype,twind,sis,hgtl_fu
 
                       call deter_sfc2(dlat_earth,dlon_earth,t4dv,idomsfc,skint,ff10,sfcr)
 
-                      write(6,*) ' radar2 ',height
                       cdata(1) = error             ! wind obs error (m/s)
                       cdata(2) = dlon              ! grid relative longitude
                       cdata(3) = dlat              ! grid relative latitude
@@ -1983,7 +1981,7 @@ subroutine read_radar(nread,ndata,nodata,infile,lunout,obstype,twind,sis,hgtl_fu
                         nodata = min(nodata+1,maxobs)  
                         usage  = zero
                         if(icuse(ikx) < 0)usage=r100
-                        if(pmot >=2 .and. usage >= 100._r_kind) rusage(ndata)=.false.
+                        if(usage >= 100._r_kind) rusage(ndata)=.false.
 
                         call deter_sfc2(dlat_earth,dlon_earth,t4dv,idomsfc,skint,ff10,sfcr)
 
@@ -2481,7 +2479,7 @@ subroutine read_radar(nread,ndata,nodata,infile,lunout,obstype,twind,sis,hgtl_fu
               if(ncnumgrp(ikx) > 0 )then                     ! cross validation on
                  if(mod(ndata,ncnumgrp(ikx))== ncgroup(ikx)-1)usage=ncmiter(ikx)
               end if
-              if(pmot >=2 .and. usage >= 100._r_kind) rusage(ndata)=.false.
+              if(usage >= 100._r_kind) rusage(ndata)=.false.
 
               call deter_zsfc_model(dlat,dlon,zsges)
 
@@ -3003,7 +3001,7 @@ subroutine read_radar(nread,ndata,nodata,infile,lunout,obstype,twind,sis,hgtl_fu
                     if(ncnumgrp(ikx) > 0 )then                     ! cross validation on
                        if(mod(ndata,ncnumgrp(ikx))== ncgroup(ikx)-1)usage=ncmiter(ikx)
                     end if
-                    if(pmot >=2 .and. usage >= 100._r_kind) rusage(ndata)=.false.
+                    if(usage >= 100._r_kind) rusage(ndata)=.false.
 
                     call deter_zsfc_model(dlat,dlon,zsges)
 

src/gsi/read_radar_wind_ascii.f90

@@ -201,7 +201,7 @@ subroutine read_radar_wind_ascii(nread,ndata,nodata,infile,lunout,obstype,sis,hg
                   clat0,slat0,dlat,dlon,thiserr,thislon,thislat, &
                   rlonloc,rlatloc,rlonglob,rlatglob,timeb,rad_per_meter
   real(r_kind) :: azm,cosazm_earth,sinazm_earth,cosazm,sinazm  
-  real(r_kind) :: radartwindow
+  real(r_kind) :: radartwindow,usage
   real(r_kind) :: rmins_an,rmins_ob                                                     
   real(r_kind),allocatable,dimension(:,:):: cdata_all
   real(r_double) rstation_id
@@ -516,7 +516,8 @@ subroutine read_radar_wind_ascii(nread,ndata,nodata,infile,lunout,obstype,sis,hg
                  save_all=.false.
                  if(pmot /= 2 .and. pmot /= 0) save_all=.true.
 
-                 if(pmot >= 2 .and. abs(icuse(ikx)) /= 1)rusage(ndata+1)=.false.
+                 usage = zero
+                 if(abs(icuse(ikx)) /= 1)usage=r100
 
                  if(ithin == 1)then
                    if(zflag == 0)then
@@ -595,6 +596,7 @@ subroutine read_radar_wind_ascii(nread,ndata,nodata,infile,lunout,obstype,sis,hg
                  cdata_all(22,iout)=two                            ! Level 2 data
 
                  if(doradaroneob .and. (cdata_all(5,iout) > -99_r_kind) ) exit volumes
+                 if(usage >= r100)rusage(iout)=.false.
 
                end do azms  !j
             else