Update to reproduce initial results (except hafs_3d)

src/enkf/letkf.f90

@@ -276,24 +276,24 @@ subroutine letkf_update()
 
 ! Update ensemble on model grid.
 ! Loop for each horizontal grid points on this task.
-! !$omp parallel do schedule(dynamic) default(none) private(npt,nob,nobsl, &
-! !$omp                  nobsl2,ngrd1,corrlength,ens_tmp,coslat, &
-! !$omp                  nf,vdist,obens,indxassim,indxob,maxdfs, &
-! !$omp                  nn,hxens,wts_ensmean,dfs,rdiag,dep,rloc,i, &
-! !$omp                  oindex,deglat,dist,corrsq,nb,nlev,nanal,sresults, &
-! !$omp                  wts_ensperts,pa,trpa,trpa_raw) shared(anal_ob, &
-! !$omp                  anal_ob_modens,anal_chunk,obsprd_post,obsprd_prior, &
-! !$omp                  oberrvar,oberrvaruse,nobsl_max,grdloc_chunk, &
-! !$omp                  obloc,corrlengthnh,corrlengthsh,corrlengthtr,&
-! !$omp                  vlocal_evecs,vlocal,oblnp,lnp_chunk,lnsigl,corrlengthsq,&
-! !$omp                  getkf,denkf,getkf_inflation,ensmean_chunk,ob,ensmean_ob, &
-! !$omp                  nproc,numptsperproc,nnmax,r_nanalsm1,kdtree_obs2,kdobs, &
-! !$omp                  mincorrlength_factsq,robs_local,coslats_local, &
-! !$omp                  lupd_obspace_serial,eps,dfs_sort,nanals,index_pres,&
-! !$omp  neigv,nlevs,lonsgrd,latsgrd,nobstot,nens,ncdim,nbackgrounds,indxproc,rad2deg) &
-! !$omp  reduction(+:t1,t2,t3,t4,t5) &
-! !$omp  reduction(max:nobslocal_max) &
-! !$omp  reduction(min:nobslocal_min) 
+!$omp parallel do schedule(dynamic) default(none) private(npt,nob,nobsl, &
+!$omp                  nobsl2,ngrd1,corrlength,ens_tmp,coslat, &
+!$omp                  nf,vdist,obens,indxassim,indxob,maxdfs, &
+!$omp                  nn,hxens,wts_ensmean,dfs,rdiag,dep,rloc,i, &
+!$omp                  oindex,deglat,dist,corrsq,nb,nlev,nanal,sresults, &
+!$omp                  wts_ensperts,pa,trpa,trpa_raw) shared(anal_ob, &
+!$omp                  anal_ob_modens,anal_chunk,obsprd_post,obsprd_prior, &
+!$omp                  oberrvar,oberrvaruse,nobsl_max,grdloc_chunk, &
+!$omp                  obloc,corrlengthnh,corrlengthsh,corrlengthtr,&
+!$omp                  vlocal_evecs,vlocal,oblnp,lnp_chunk,lnsigl,corrlengthsq,&
+!$omp                  getkf,denkf,getkf_inflation,ensmean_chunk,ob,ensmean_ob, &
+!$omp                  nproc,numptsperproc,nnmax,r_nanalsm1,kdtree_obs2,kdobs, &
+!$omp                  mincorrlength_factsq,robs_local,coslats_local, &
+!$omp                  lupd_obspace_serial,eps,dfs_sort,nanals,index_pres,&
+!$omp  neigv,nlevs,lonsgrd,latsgrd,nobstot,nens,ncdim,nbackgrounds,indxproc,rad2deg) &
+!$omp  reduction(+:t1,t2,t3,t4,t5) &
+!$omp  reduction(max:nobslocal_max) &
+!$omp  reduction(min:nobslocal_min) 
 grdloop: do npt=1,numptsperproc(nproc+1)
 
    t1 = mpi_wtime()

src/gsi/gsdcloudanalysis.F90

@@ -64,8 +64,8 @@ subroutine  gsdcloudanalysis(mype)
 !_____________________________________________________________________
 !
 ! 
-  use constants, only: zero,one
-  use constants, only: h1000
+  use constants, only: zero,one,rad2deg,fv
+  use constants, only: rd_over_cp,h1000
   use kinds,   only: r_single,i_kind, r_kind
   use gridmod, only: pt_ll,eta1_ll,aeta1_ll,eta2_ll,aeta2_ll
   use gridmod, only: regional,wrf_mass_regional,regional_time

src/gsi/intrw.f90

@@ -127,23 +127,17 @@ subroutine intrw_(rwhead,rval,sval)
   ier=0
   call gsi_bundlegetpointer(sval,'u',su,istatus);ier=istatus+ier
   call gsi_bundlegetpointer(sval,'v',sv,istatus);ier=istatus+ier
-  call gsi_bundlegetpointer(sval,'w',sw,istatus)
-  if (if_use_w_vr.and.istatus==0) then
-     include_w=.true.
-  else
-     include_w=.false.
-  end if
   call gsi_bundlegetpointer(rval,'u',ru,istatus);ier=istatus+ier
   call gsi_bundlegetpointer(rval,'v',rv,istatus);ier=istatus+ier
-  call gsi_bundlegetpointer(rval,'w',rw,istatus)
-  if (if_use_w_vr.and.istatus==0) then
-     include_w=.true.
-  else
-     include_w=.false.
-  end if
 
   if(ier/=0)return
 
+  include_w=.false.
+  call gsi_bundlegetpointer(sval,'w',sw,istatus)
+  if (if_use_w_vr.and.istatus==0) then
+     call gsi_bundlegetpointer(rval,'w',rw,istatus)
+     if(istatus == 0)include_w=.true.
+  end if
 
   !rwptr => rwhead
   rwptr => rwNode_typecast(rwhead)

src/gsi/read_prepbufr.f90

@@ -3010,6 +3010,16 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
   call closbf(lunin)
   close(lunin)
 
+! Apply hilbert curve for cross validation if requested
+
+  if(lhilbert) then
+       call apply_hilbertcurve(ndata,obstype,cdata_all(thisobtype_usage,1:ndata))   
+
+     do i=1,ndata
+        if(cdata_all(thisobtype_usage,i) >= 100._r_kind) rusage(i) = .false.
+     end do
+  end if
+
   nxdata=ndata
   ndata=0
   if(nxdata > 0)then
@@ -3057,10 +3067,6 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
   end if
   deallocate(rusage,rthin)
 
-! Apply hilbert curve for cross validation if requested
-
-  if(lhilbert) &
-       call apply_hilbertcurve(ndata,obstype,cdata_all(thisobtype_usage,1:ndata))   
 
 ! the following is gettin the types which will be applied hilbert curve to
 !  estimate the density

src/gsi/read_radar.f90

@@ -99,14 +99,14 @@ subroutine read_radar(nread,ndata,nodata,infile,lunout,obstype,twind,sis,hgtl_fu
   use convinfo, only: nconvtype,ctwind, &
       ncmiter,ncgroup,ncnumgrp,icuse,ictype,ioctype,ithin_conv,rmesh_conv,pmesh_conv,pmot_conv
   use convthin, only: make3grids,map3grids_m,del3grids,use_all
-  use deter_sfc_mod, only: deter_sfc2,deter_zsfc_model
   use mpimod, only: npe
   use gsi_io, only: verbose
   use mpimod, only: mype
   use directDA_radaruse_mod, only: oe_rw, lvldbg, refl_lowbnd_rw
   use directDA_radaruse_mod, only: l_correct_azmu, l_correct_tilt, i_correct_tilt, &
                                l_azm_east1st, l_plt_diag_rw
   use directDA_radaruse_mod, only: l_use_rw_columntilt
+  use deter_sfc_mod, only: deter_sfc2,deter_zsfc_model
 
   implicit none 
 
@@ -355,8 +355,6 @@ subroutine read_radar(nread,ndata,nodata,infile,lunout,obstype,twind,sis,hgtl_fu
   dlonmax=-huge(dlonmax)
   dlatmin=huge(dlatmin)
   dlonmin=huge(dlonmin)
-  toff=zero                     !uncertainty
-!  uncertainty means that there is an issue with defining the variable.
 
   if(ianldate > 2016092000)then
      hdrstr(2)='PTID YEAR MNTH DAYS HOUR MINU SECO CLAT CLON FLVLST ANAZ ANEL'
@@ -1299,7 +1297,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_sfc2(dlat_earth,dlon_earth,t4dv,idomsfc,skint,ff10,sfcr)
 
@@ -2178,7 +2176,6 @@ subroutine read_radar(nread,ndata,nodata,infile,lunout,obstype,twind,sis,hgtl_fu
   ibadstaheight=0
   notgood=0
   notgood0=0
-  nread=0
   ntdrvr_in=0
   ntdrvr_kept=0
   ntdrvr_thin1=0
@@ -2804,8 +2801,9 @@ subroutine read_radar(nread,ndata,nodata,infile,lunout,obstype,twind,sis,hgtl_fu
               nread=nread+1
 !             Select data every 3 km along each beam
               if(MOD(INT(tdr_obs(1,k)-tdr_obs(1,1)),3000) < 100)then
-                 if(tdr_obs(3,k) >= 800.) nmissing=nmissing+1     !xx
-                 if(tdr_obs(3,k) < 800.) then
+                 if(tdr_obs(3,k) >= 800.) then
+                    nmissing=nmissing+1     !xx
+                 else
                     ii=ii+1
                     dopbin(ii)=tdr_obs(3,k)
                     thisrange=tdr_obs(1,k)
@@ -2922,6 +2920,7 @@ subroutine read_radar(nread,ndata,nodata,infile,lunout,obstype,twind,sis,hgtl_fu
                  good=.true.
                  if(.not.good0) then
                     notgood0=notgood0+1
+                    good=.false.
                     cycle
                  end if
 !           if data is good, load into output array
@@ -3425,7 +3424,6 @@ subroutine read_radar_l2rw_novadqc(ndata,nodata,lunout,obstype,sis,nobs)
   notgood0=0
   nsuper2_in=0
   nsuper2_kept=0
-  toff=zero          !uncertainty
 
   if(loop==0) outmessage='level 2 superobs:'
 
@@ -3479,7 +3477,6 @@ subroutine read_radar_l2rw_novadqc(ndata,nodata,lunout,obstype,sis,nobs)
         timeo=thistime
         if(abs(timeo)>half ) cycle
      endif
-     t4dv = thistime                 ! uncertainty
 
 !    Get observation (lon,lat).  Compute distance from radar.
      dlat_earth=thislat
@@ -3569,12 +3566,12 @@ subroutine read_radar_l2rw_novadqc(ndata,nodata,lunout,obstype,sis,nobs)
      if(thiserr>r6 .or. thiserr<=zero) then
         ibaderror=ibaderror+1; good0=.false.
      end if
+     good=.true.
      if(.not.good0) then
         good=.false.
         notgood0=notgood0+1
         cycle
      end if
-     good=.true.
 
 !    If data is good, load into output array
      if(good) then
@@ -3588,7 +3585,6 @@ subroutine read_radar_l2rw_novadqc(ndata,nodata,lunout,obstype,sis,nobs)
         end if
         if(usage >= 100._r_kind)rusage(ndata)=.true.
 
-        call deter_zsfc_model(dlat,dlon,zsges)
         call deter_sfc2(dlat_earth,dlon_earth,t4dv,idomsfc,skint,ff10,sfcr)
 
         cdata(1) = error             ! wind obs error (m/s)
@@ -3747,7 +3743,7 @@ subroutine read_radar_l2rw(ndata,nodata,lunout,obstype,sis,nobs,hgtl_full)
   integer(i_kind) iret,kx0
   integer(i_kind) nreal,nchanl,ilat,ilon,ikx
   integer(i_kind) idomsfc
-  real(r_kind) usage,ff10,sfcr,skint,t4dvo,toff
+  real(r_kind) usage,ff10,sfcr,skint,t4dv,t4dvo,toff
   real(r_kind) eradkm,dlat_earth,dlon_earth
   real(r_kind) dlat,dlon,staheight,tiltangle,clon,slon,clat,slat
   real(r_kind) timeo,clonh,slonh,clath,slath,cdist,dist
@@ -3788,7 +3784,7 @@ subroutine read_radar_l2rw(ndata,nodata,lunout,obstype,sis,nobs,hgtl_full)
   character(4),allocatable,dimension(:):: rsite 
   integer(i_kind),allocatable,dimension(:):: ruse
   character(8) chdr2,subset 
-  real(r_double) rdisttest(n_gates_max),hdr(10),hdr2(12),rwnd0(3,n_gates_max) 
+  real(r_double) rdisttest(n_gates_max),hdr(3),hdr2(12),rwnd0(3,n_gates_max) 
   character(4) stn_id 
   equivalence (chdr2,hdr2(1)) 
   real(r_kind) stn_lat,stn_lon,stn_hgt,stn_az,stn_el,t,range,vrmax,vrmin,aactual,a43,b,c,selev0,celev0,thistiltr,epsh,h,ha,rlonloc,rlatloc
@@ -3797,7 +3793,7 @@ subroutine read_radar_l2rw(ndata,nodata,lunout,obstype,sis,nobs,hgtl_full)
   real(r_kind):: relm,srlm,crlm,sph,cph,cc,anum,denom 
   real(r_kind) :: rmesh,xmesh,zmesh,dx,dy,dx1,dy1,w00,w01,w10,w11 
   real(r_kind), allocatable, dimension(:) :: zl_thin 
-  integer(i_kind) :: ithin,zflag,nlevz,icntpnt,klon1,klat1,kk,klatp1,klonp1 
+  integer(i_kind) :: ithin,zflag,nlevz,klon1,klat1,kk,klatp1,klonp1 
   real(r_kind),dimension(nsig):: hges,zges 
   real(r_kind) sin2,termg,termr,termrg,zobs 
   integer(i_kind) iout,ntdrvr_thin2 
@@ -3854,7 +3850,6 @@ subroutine read_radar_l2rw(ndata,nodata,lunout,obstype,sis,nobs,hgtl_full)
   nmrecs=0
   irec=0
   errzmax=zero
-  toff=zero                !uncertainty
 
   timemax=-huge(timemax)
   timemin=huge(timemin)
@@ -3877,7 +3872,6 @@ subroutine read_radar_l2rw(ndata,nodata,lunout,obstype,sis,nobs,hgtl_full)
   nsuper2_in=0
   nsuper2_kept=0
   ntdrvr_thin2=0
-  icntpnt=0
   nout=0
   if(loop==0) outmessage='level 2 superobs:'
   rmesh=radar_rmesh
@@ -3965,11 +3959,10 @@ subroutine read_radar_l2rw(ndata,nodata,lunout,obstype,sis,nobs,hgtl_full)
       stn_lat=hdr2(2)
       stn_lon=hdr2(3)
       stn_hgt=hdr2(4)+hdr2(5)
-      call ufbint(inbufr,hdr,10,1,levs, &
-           'SSTN YEAR MNTH DAYS HOUR MINU SECO ANAZ ANEL QCRW')
+      call ufbint(inbufr,hdr,3,1,levs,'ANAZ ANEL QCRW')
       nradials_in=nradials_in+1
-      stn_az=r90-hdr(8)
-      stn_el=hdr(9)
+      stn_az=r90-hdr(1)
+      stn_el=hdr(2)
       call ufbint(inbufr,rwnd0,3,n_gates_max,n_gates,'DIST125M DMVR DVSW')
       do i=1,n_gates
         range=distfact*rwnd0(1,i)
@@ -4238,7 +4231,9 @@ subroutine read_radar_l2rw(ndata,nodata,lunout,obstype,sis,nobs,hgtl_full)
             usage=r100
           end if
 
-          call deter_zsfc_model(dlat,dlon,zsges)
+!    Get information from surface file necessary for conventional data here
+!         call deter_zsfc_model(dlat,dlon,zsges)
+!         call deter_sfc2(dlat_earth,dlon_earth,t4dv,idomsfc,skint,ff10,sfcr)
 
           nsuper2_kept=nsuper2_kept+1
           cdata(1) = error             ! wind obs error (m/s)
@@ -4328,6 +4323,7 @@ subroutine read_radar_l2rw(ndata,nodata,lunout,obstype,sis,nobs,hgtl_full)
   write(6,*)'READ_RADAR_L2RW: dlatmin,max,dlonmin,max=',dlatmin,dlatmax,dlonmin,dlonmax
 
 ! Write observation to scratch file
+  deallocate(rusage,rthin)
   call count_obs(ndata,maxdat,ilat,ilon,cdata_all,nobs)
   write(lunout) obstype,sis,nreal,nchanl,ilat,ilon
   write(lunout) ((cdata_all(k,i),k=1,maxdat),i=1,ndata)

src/gsi/read_satwnd.f90

@@ -706,7 +706,6 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis
                  enddo
               endif
            endif
-           write(6,*) ioctype(nc),ictype(nc),rmesh,pflag,nlevp,pmesh,nc,pmot,ptime
            write(6,'(a52,a16,I5,f10.2,2i5,f10.2,i5,i5,f10.2)') &
                    ' READ_SATWND: ictype(nc),rmesh,pflag,nlevp,pmesh,nc ', &
                    ioctype(nc),ictype(nc),rmesh,pflag,nlevp,pmesh,nc,pmot,ptime