From f76d872859257bfa614d6f6395132700ca6b1b77 Mon Sep 17 00:00:00 2001
From: jderber-NOAA <75998838+jderber-NOAA@users.noreply.github.com>
Date: Tue, 17 Oct 2023 10:43:33 -0400
Subject: [PATCH] Replacement of float command with real command (#631)

---
 src/gsi/aniso_ens_util.f90                 |  44 +-
 src/gsi/anisofilter.f90                    |  50 +--
 src/gsi/anisofilter_glb.f90                |  26 +-
 src/gsi/atms_spatial_average_mod.f90       |   2 +-
 src/gsi/balmod.f90                         |  10 +-
 src/gsi/berror.f90                         |   4 +-
 src/gsi/bkgvar_rewgt.f90                   |   2 +-
 src/gsi/buddycheck_mod.f90                 |   2 +-
 src/gsi/calc_fov_crosstrk.f90              |   2 +-
 src/gsi/compact_diffs.f90                  |  54 +--
 src/gsi/compute_derived.f90                |   2 +-
 src/gsi/compute_qvar3d.f90                 |   4 +-
 src/gsi/cplr_get_fv3_regional_ensperts.f90 |   4 +-
 src/gsi/cplr_get_pseudo_ensperts.f90       |  10 +-
 src/gsi/cplr_get_wrf_mass_ensperts.f90     |   4 +-
 src/gsi/cplr_get_wrf_nmm_ensperts.f90      |   2 +-
 src/gsi/cplr_gfs_ensmod.f90                |  13 +-
 src/gsi/cplr_read_wrf_nmm_guess.f90        |   8 +-
 src/gsi/cplr_wrf_netcdf_interface.f90      |   4 +-
 src/gsi/deter_sfc_mod.f90                  |  22 +-
 src/gsi/ens_spread_mod.f90                 |   2 +-
 src/gsi/general_commvars_mod.f90           |  20 +-
 src/gsi/general_read_fv3atm.f90            |   3 +-
 src/gsi/general_read_gfsatm.f90            |  33 +-
 src/gsi/general_read_nemsaero.f90          |   3 +-
 src/gsi/general_spectral_transforms.f90    |  13 +-
 src/gsi/general_sub2grid_mod.f90           |   4 +-
 src/gsi/gengrid_vars.f90                   |   4 +-
 src/gsi/gesinfo.F90                        |  11 +-
 src/gsi/get_gefs_ensperts_dualres.f90      |   4 +-
 src/gsi/get_gefs_for_regional.f90          |  10 +-
 src/gsi/get_nmmb_ensperts.f90              |   2 +-
 src/gsi/gfs_stratosphere.f90               |  12 +-
 src/gsi/grdcrd.f90                         |   4 +-
 src/gsi/gsd_update_mod.f90                 |   2 +-
 src/gsi/gsdcloudlib_pseudoq_mod.f90        |   6 +-
 src/gsi/guess_grids.F90                    |   6 +-
 src/gsi/hdraobmod.f90                      |   2 +-
 src/gsi/hybrid_ensemble_isotropic.F90      |   2 +-
 src/gsi/intjcmod.f90                       |   4 +-
 src/gsi/intrp2a.f90                        |  12 +-
 src/gsi/m_berror_stats.f90                 |   4 +-
 src/gsi/m_extOzone.F90                     |   4 +-
 src/gsi/mod_fv3_lola.f90                   |   4 +-
 src/gsi/mod_vtrans.f90                     |   2 +-
 src/gsi/mod_wrfmass_to_a.f90               |   6 +-
 src/gsi/mp_compact_diffs_mod1.f90          |  24 +-
 src/gsi/ncepgfs_ghg.f90                    |   4 +-
 src/gsi/ncepgfs_io.f90                     |  10 +-
 src/gsi/ncepnems_io.f90                    |  19 +-
 src/gsi/netcdfgfs_io.f90                   |   4 +-
 src/gsi/nlmsas_ad.f90                      |   6 +-
 src/gsi/polcarf.f90                        |   4 +-
 src/gsi/prewgt.f90                         |   4 +-
 src/gsi/prewgt_reg.f90                     |  12 +-
 src/gsi/q_diag.f90                         |   2 +-
 src/gsi/radinfo.f90                        |   6 +-
 src/gsi/raflib.f90                         |   4 +-
 src/gsi/rdgrbsst.f90                       |   8 +-
 src/gsi/read_airs.f90                      |   2 +-
 src/gsi/read_amsr2.f90                     |   2 +-
 src/gsi/read_atms.f90                      |   8 +-
 src/gsi/read_bufrtovs.f90                  |   6 +-
 src/gsi/read_cris.f90                      |   8 +-
 src/gsi/read_dbz_netcdf.f90                | 460 ++++++++++-----------
 src/gsi/read_files.f90                     |  19 +-
 src/gsi/read_gfs_ozone_for_regional.f90    |   4 +-
 src/gsi/read_gmi.f90                       |   6 +-
 src/gsi/read_goesglm.f90                   |  34 +-
 src/gsi/read_goesndr.f90                   |   2 +-
 src/gsi/read_iasi.f90                      |   6 +-
 src/gsi/read_lidar.f90                     |   2 +-
 src/gsi/read_nsstbufr.f90                  |   2 +-
 src/gsi/read_ozone.f90                     |   4 +-
 src/gsi/read_pblh.f90                      |   2 +-
 src/gsi/read_prepbufr.f90                  |   2 +-
 src/gsi/read_radar_wind_ascii.f90          | 242 +++++------
 src/gsi/read_saphir.f90                    |  10 +-
 src/gsi/read_wcpbufr.f90                   |   2 +-
 src/gsi/reorg_metar_cloud.f90              |  10 +-
 src/gsi/rfdpar.f90                         |   2 +-
 src/gsi/satthin.F90                        |   4 +-
 src/gsi/setupbend.f90                      |   6 +-
 src/gsi/setupdw.f90                        |   4 +-
 src/gsi/setuplag.f90                       |   2 +-
 src/gsi/setuplight.f90                     |   2 +-
 src/gsi/setupoz.f90                        |   6 +-
 src/gsi/setuppcp.f90                       |   4 +-
 src/gsi/setupq.f90                         |   2 +-
 src/gsi/setupref.f90                       |   6 +-
 src/gsi/setupspd.f90                       |   2 +-
 src/gsi/setupt.f90                         |   6 +-
 src/gsi/setuptcp.f90                       |   4 +-
 src/gsi/setupw.f90                         |   4 +-
 src/gsi/setupwspd10m.f90                   |   2 +-
 src/gsi/sfcobsqc.f90                       |   2 +-
 src/gsi/smoothzrf.f90                      |   2 +-
 src/gsi/ssmis_spatial_average_mod.f90      |   2 +-
 src/gsi/statsco.f90                        |   8 +-
 src/gsi/statsconv.f90                      |  74 ++--
 src/gsi/statsoz.f90                        |   4 +-
 src/gsi/statspcp.f90                       |   2 +-
 src/gsi/statsrad.f90                       |   2 +-
 src/gsi/stpjcmod.f90                       |   2 +-
 src/gsi/support_2dvar.f90                  |  22 +-
 src/gsi/tcv_mod.f90                        |   8 +-
 src/gsi/tintrp2a.f90                       |  16 +-
 src/gsi/tintrp3.f90                        |  12 +-
 src/gsi/wind_fft.f90                       |   8 +-
 src/gsi/write_fv3_spread.f90               |   6 +-
 src/gsi/write_incr.f90                     |   6 +-
 111 files changed, 817 insertions(+), 803 deletions(-)

diff --git a/src/gsi/aniso_ens_util.f90 b/src/gsi/aniso_ens_util.f90
index f118bee40f..43f81216e4 100644
--- a/src/gsi/aniso_ens_util.f90
+++ b/src/gsi/aniso_ens_util.f90
@@ -122,8 +122,8 @@ subroutine ens_uv_to_psichi(u,v,truewind)
         do j=1,nlon
            rlon=region_lon(i,j)
            rlat=region_lat(i,j)
-           dlon=float(j)*one
-           dlat=float(i)*one
+           dlon=real(j,r_kind)
+           dlat=real(i,r_kind)
            ue=u(i,j)
            ve=v(i,j)
            call rotate_wind_ll2xy(ue,ve,ug,vg,rlon,dlon,dlat)
@@ -440,13 +440,13 @@ subroutine ens_intpcoeffs_reg(ngrds,igbox,iref,jref,igbox0f,ensmask,enscoeff,gbl
               yg=rlat+90._r_kind+one
            end if
 
-           dxg=xg-float(floor(xg))
-           dyg=yg-float(floor(yg))
+           dxg=xg-real(floor(xg),r_kind)
+           dyg=yg-real(floor(yg),r_kind)
            dxg1=one-dxg
            dyg1=one-dyg
 
-           if (xg>=one .and. xg<=float(jxp) .and.  &
-               yg>=one .and. yg<=float(iy) ) then
+           if (xg>=one .and. xg<=real(jxp,r_kind) .and.  &
+               yg>=one .and. yg<=real(iy,r_kind) ) then
 
               enscoeff(1,i,j,kg)=dxg1*dyg1
               enscoeff(2,i,j,kg)=dxg1*dyg
@@ -479,9 +479,9 @@ subroutine ens_intpcoeffs_reg(ngrds,igbox,iref,jref,igbox0f,ensmask,enscoeff,gbl
      endif
 
      do j=1,iy
-        yg=float(j)*one
+        yg=real(j,r_kind)
         do i=1,jx
-           xg=float(i)*one
+           xg=real(i,r_kind)
            call w3fb12(xg,yg,alat1,elon1,ds,elonv,alatan,rlat,rlon,ierr8)
            rlon=rlon/rad2deg
            rlat=rlat/rad2deg
@@ -620,34 +620,34 @@ subroutine ens_intpcoeffs_reg(ngrds,igbox,iref,jref,igbox0f,ensmask,enscoeff,gbl
      igbox(2,kg)=iimax0(kg)
      igbox(3,kg)=jjmin0(kg)
      igbox(4,kg)=jjmax0(kg)
-     igbox0f(1,kg)=one+float((igbox(1,kg)-1))/pf2aP1%grid_ratio_lat + ijadjust
-     igbox0f(2,kg)=one+float((igbox(2,kg)-1))/pf2aP1%grid_ratio_lat - ijadjust
-     igbox0f(3,kg)=one+float((igbox(3,kg)-1))/pf2aP1%grid_ratio_lon + ijadjust
-     igbox0f(4,kg)=one+float((igbox(4,kg)-1))/pf2aP1%grid_ratio_lon - ijadjust
+     igbox0f(1,kg)=one+real(igbox(1,kg)-1,r_kind)/pf2aP1%grid_ratio_lat + ijadjust
+     igbox0f(2,kg)=one+real(igbox(2,kg)-1,r_kind)/pf2aP1%grid_ratio_lat - ijadjust
+     igbox0f(3,kg)=one+real(igbox(3,kg)-1,r_kind)/pf2aP1%grid_ratio_lon + ijadjust
+     igbox0f(4,kg)=one+real(igbox(4,kg)-1,r_kind)/pf2aP1%grid_ratio_lon - ijadjust
   end do
 
 !==> compute blending functions
 
   do i=1,pf2aP1%nlatf
-     dist1=float(igbox0f(1,1)-i)
-     dist2=float(i-igbox0f(2,1))
+     dist1=real(igbox0f(1,1)-i,r_kind)
+     dist2=real(i-igbox0f(2,1),r_kind)
      gblend_b(i,1)=half*(one-tanh(dist1)) !relax to zero
      gblend_t(i,1)=half*(one-tanh(dist2)) !outside 212 grid
 
-     dist1=float(igbox0f(1,2)-i)
-     dist2=float(i-igbox0f(2,2))
+     dist1=real(igbox0f(1,2)-i,r_kind)
+     dist2=real(i-igbox0f(2,2),r_kind)
      gblend_b(i,2)=half*(one-tanh(dist1)) !relax to zero
      gblend_t(i,2)=half*(one-tanh(dist2)) !outside 221 grid
   end do
 
   do j=1,pf2aP1%nlonf
-     dist1=float(igbox0f(3,1)-j)
-     dist2=float(j-igbox0f(4,1))
+     dist1=real(igbox0f(3,1)-j,r_kind)
+     dist2=real(j-igbox0f(4,1),r_kind)
      gblend_l(j,1)=half*(one-tanh(dist1)) !relax to zero
      gblend_r(j,1)=half*(one-tanh(dist2)) !outside 212 grid
  
-     dist1=float(igbox0f(3,2)-j)
-     dist2=float(j-igbox0f(4,2))
+     dist1=real(igbox0f(3,2)-j,r_kind)
+     dist2=real(j-igbox0f(4,2),r_kind)
      gblend_l(j,2)=half*(one-tanh(dist1)) !relax to zero
      gblend_r(j,2)=half*(one-tanh(dist2)) !outside 221 grid
   end do
@@ -1141,10 +1141,10 @@ subroutine ens_fill(ur,na,nb,u,nxx,ny,itap,no_wgt_in)
   no_wgt=.false.
   if(no_wgt_in) no_wgt=.true.
 
-  pionp1=four*atan(one)/float(itap+1)
+  pionp1=four*atan(one)/real(itap+1,r_kind)
 
   do i=1,itap
-     xi=float(i)
+     xi=real(i,r_kind)
      wt(i)=half+half*cos(pionp1*xi)
   enddo
 
diff --git a/src/gsi/anisofilter.f90 b/src/gsi/anisofilter.f90
index ec05d191ba..c05c764a05 100755
--- a/src/gsi/anisofilter.f90
+++ b/src/gsi/anisofilter.f90
@@ -596,7 +596,7 @@ subroutine anprewgt_reg(mype)
         do i=indices%ips,indices%ipe
            l =max(min(int(rllatf(i,j)),mlat),1)
            lp=min((l+1),mlat)
-           dl2=rllatf(i,j)-float(l)
+           dl2=rllatf(i,j)-real(l,r_kind)
            dl1=one-dl2
            if (ivar <= nrf) then
               if (nrf_3d(ivar)) then
@@ -1056,7 +1056,7 @@ subroutine get_aspect_reg_pt(mype)
               asp3=scalex3*asp3
            endif
 
-           rk1=float(k1-44)
+           rk1=real(k1-44,r_kind)
            fblend=half*(one-tanh(rk1))! one
 
            if (nvar_id(k) /= nrf3_loc(nrf3_q)) then
@@ -1126,7 +1126,7 @@ subroutine fact_qopt2(factk,rh,kvar)
   d  =20.0_r_kind * rh + one
   n  =int(d)
   np =n+1
-  dn2=d-float(n)
+  dn2=d-real(n,r_kind)
   dn1=one-dn2
   n =min0(max(1,n) ,25)
   np=min0(max(1,np),25)
@@ -2407,7 +2407,7 @@ subroutine read_bckgstats(mype)
      do k=1,nsig
         vzimax(k,n)=maxval(one/vz(k,0:mlat+1,n))
         vzimin(k,n)=minval(one/vz(k,0:mlat+1,n))
-        vziavg(k,n)=sum((one/vz(k,0:mlat+1,n)))/float(mlat+2)
+        vziavg(k,n)=sum((one/vz(k,0:mlat+1,n)))/real(mlat+2,r_kind)
      end do
      if(print_verbose) then
         do k=1,nsig
@@ -2428,13 +2428,13 @@ subroutine read_bckgstats(mype)
 
      do n=1,nrf3
         do k=1,nsig
-           corzavg(k,n)=sum(corz(1:mlat,k,n))/float(mlat)
-           hwllavg(k,n)=sum(hwll(0:mlat+1,k,n))/float(mlat+2)
+           corzavg(k,n)=sum(corz(1:mlat,k,n))/real(mlat,r_kind)
+           hwllavg(k,n)=sum(hwll(0:mlat+1,k,n))/real(mlat+2,r_kind)
         end do
      end do
      do n=1,nvars-nrf3
-        corpavg(n)=sum(corp(1:mlat,n))/float(mlat)
-        hwllpavg(n)=sum(hwllp(0:mlat+1,n))/float(mlat+2)
+        corpavg(n)=sum(corp(1:mlat,n))/real(mlat,r_kind)
+        hwllpavg(n)=sum(hwllp(0:mlat+1,n))/real(mlat+2,r_kind)
      end do
 
      do j=1,mlat
@@ -2869,7 +2869,7 @@ subroutine isotropic_scales(scale1,scale2,scale3,k)
                     else
                        l =max(min(int(rllatf(i,j)),mlat),1)
                        lp=min((l+1),mlat)
-                       dl2=rllatf(i,j)-float(l)
+                       dl2=rllatf(i,j)-real(l,r_kind)
                        dl1=one-dl2
                        hwll_loc=dl1*hwll(l,k1,n)+dl2*hwll(lp,k1,n)
                     end if
@@ -2886,7 +2886,7 @@ subroutine isotropic_scales(scale1,scale2,scale3,k)
 
                  l =max(min(int(rllatf(i,j)),mlat),1)
                  lp=min((l+1),mlat)
-                 dl2=rllatf(i,j)-float(l)
+                 dl2=rllatf(i,j)-real(l,r_kind)
                  dl1=one-dl2
                  hwll_loc=cc*(dl1*hwllp(l,n)+dl2*hwllp(lp,n))
                  scale3(i,j)=one
@@ -2903,7 +2903,7 @@ subroutine isotropic_scales(scale1,scale2,scale3,k)
 
                  l =max(min(int(rllatf(i,j)),mlat),1)
                  lp=min((l+1),mlat)
-                 dl2=rllatf(i,j)-float(l)
+                 dl2=rllatf(i,j)-real(l,r_kind)
                  dl1=one-dl2
                  hwll_loc=cc*(dl1*hwllp(l,nn)+dl2*hwllp(lp,nn))
                  scale3(i,j)=one
@@ -3027,7 +3027,7 @@ subroutine get_theta_corrl_lenghts(mype)
      mcount0=lon2*lat2! It's OK to count buffer points
      call mpi_allreduce(pbar4a,pbar4(k),1,mpi_real8,mpi_sum,mpi_comm_world,ierror)
      call mpi_allreduce(mcount0,mcount,1,mpi_integer4,mpi_sum,mpi_comm_world,ierror)
-     pbar4(k)=pbar4(k)/float(mcount)
+     pbar4(k)=pbar4(k)/real(mcount,r_kind)
      if(print_verbose) write(6,*)'in get_theta_corrl_lenghts,k,pbar4=',k,pbar4(k)
      call w3fa03(pbar4(k),hgt4(k),tbar4(k),thetabar4(k))
   end do
@@ -3881,15 +3881,15 @@ subroutine get_aspect_reg_ens(mype)
  
         do j=1,pf2aP1%nlonf
            do i=1,pf2aP1%nlatf
-              ensv(i,j,k,1)=(ensv(i,j,k,1)+ensv(i,j,k,2)+ensv(i,j,k,3))/sqrt(float(nt1))
-              ensv(i,j,k,2)=              (ensv(i,j,k,2)+ensv(i,j,k,3))/sqrt(float(nt2))
-              ensv(i,j,k,3)=                             ensv(i,j,k,3) /sqrt(float(nt3))
+              ensv(i,j,k,1)=(ensv(i,j,k,1)+ensv(i,j,k,2)+ensv(i,j,k,3))/sqrt(real(nt1,r_kind))
+              ensv(i,j,k,2)=              (ensv(i,j,k,2)+ensv(i,j,k,3))/sqrt(real(nt2,r_kind))
+              ensv(i,j,k,3)=                             ensv(i,j,k,3) /sqrt(real(nt3,r_kind))
 
               if( ibldani==0 .or. ibldani==2 .or. ibldani==3 ) then
                  do m=1,6
-                    c(m,1)=(aniasp(m,i,j,k,1)+aniasp(m,i,j,k,2)+aniasp(m,i,j,k,3))/float(nt1)
-                    c(m,2)=                  (aniasp(m,i,j,k,2)+aniasp(m,i,j,k,3))/float(nt2)
-                    c(m,3)=                                     aniasp(m,i,j,k,3) /float(nt3)
+                    c(m,1)=(aniasp(m,i,j,k,1)+aniasp(m,i,j,k,2)+aniasp(m,i,j,k,3))/real(nt1,r_kind)
+                    c(m,2)=                  (aniasp(m,i,j,k,2)+aniasp(m,i,j,k,3))/real(nt2,r_kind)
+                    c(m,3)=                                     aniasp(m,i,j,k,3) /real(nt3,r_kind)
                  end do
                  do igd=1,3
                     qlx=max(qlxmin(ivar,k1),ensv(i,j,k,igd))
@@ -3906,9 +3906,9 @@ subroutine get_aspect_reg_ens(mype)
                  end do
               else if(ibldani==1) then
                  do m=1,6
-                    aniasp(m,i,j,k,1)=(aniasp(m,i,j,k,1)+aniasp(m,i,j,k,2)+aniasp(m,i,j,k,3))/float(nt1)
-                    aniasp(m,i,j,k,2)=                  (aniasp(m,i,j,k,2)+aniasp(m,i,j,k,3))/float(nt2)
-                    aniasp(m,i,j,k,3)=                                     aniasp(m,i,j,k,3) /float(nt3)
+                    aniasp(m,i,j,k,1)=(aniasp(m,i,j,k,1)+aniasp(m,i,j,k,2)+aniasp(m,i,j,k,3))/real(nt1,r_kind)
+                    aniasp(m,i,j,k,2)=                  (aniasp(m,i,j,k,2)+aniasp(m,i,j,k,3))/real(nt2,r_kind)
+                    aniasp(m,i,j,k,3)=                                     aniasp(m,i,j,k,3) /real(nt3,r_kind)
                  end do
                  smax=real(maxval(ensv(i,j,k,1:3)),r_kind)
                  aensv(1,k)=aensv(1,k)+max(smax                  ,qlxmin(ivar,k1))/nlatlonf
@@ -5326,7 +5326,7 @@ subroutine get2berr_reg_subdomain_option(mype)
 
            l=max(min(int(rllatf(i,j)),mlat),1)
            lp=min((l+1),mlat)
-           dl2=rllatf(i,j)-float(l)
+           dl2=rllatf(i,j)-real(l,r_kind)
            dl1=one-dl2
            if (ivar <= nrf) then
               if (nrf_3d(ivar)) then
@@ -6520,7 +6520,7 @@ subroutine isotropic_scales_subdomain_option(scale1,scale2,scale3,k,mype)
                   else
                     l=int(rllat(iglob,jglob))
                     lp=l+1
-                    dl2=rllat(iglob,jglob)-float(l)
+                    dl2=rllat(iglob,jglob)-real(l,r_kind)
                     dl1=one-dl2
                     hwll_loc=dl1*hwll(l,k1,n)+dl2*hwll(lp,k1,n)
                     scale3(i,j)=one/vz(k1,l,n)
@@ -6536,7 +6536,7 @@ subroutine isotropic_scales_subdomain_option(scale1,scale2,scale3,k,mype)
 
                  l=int(rllat(iglob,jglob))
                  lp=l+1
-                 dl2=rllat(iglob,jglob)-float(l)
+                 dl2=rllat(iglob,jglob)-real(l,r_kind)
                  dl1=one-dl2
                  hwll_loc=cc*(dl1*hwllp(l,n)+dl2*hwllp(lp,n))
                  scale3(i,j)=one
@@ -6553,7 +6553,7 @@ subroutine isotropic_scales_subdomain_option(scale1,scale2,scale3,k,mype)
 
                  l=int(rllat(iglob,jglob))
                  lp=l+1
-                 dl2=rllat(iglob,jglob)-float(l)
+                 dl2=rllat(iglob,jglob)-real(l,r_kind)
                  dl1=one-dl2
                  hwll_loc=cc*(dl1*hwllp(l,nn)+dl2*hwllp(lp,nn))
                  scale3(i,j)=one
diff --git a/src/gsi/anisofilter_glb.f90 b/src/gsi/anisofilter_glb.f90
index 43e67a3baa..f79b26ab79 100644
--- a/src/gsi/anisofilter_glb.f90
+++ b/src/gsi/anisofilter_glb.f90
@@ -609,7 +609,7 @@ subroutine get_stat_factk(platf,ivar,kvar,factk,rh,dvsst)
 
   l =int(platf)
   lp=l+1
-  dl2=platf-float(l)
+  dl2=platf-real(l,r_kind)
   dl1=one-dl2
   l = min(max(1,l ),mlat)
   lp= min(max(1,lp),mlat)
@@ -971,7 +971,7 @@ subroutine read_bckgstats_glb(mype)
      mcount0=lon2*lat2! It's OK to count buffer points
      call mpi_allreduce(pbar4a,pbar4(k),1,mpi_real8,mpi_sum,mpi_comm_world,ierror)
      call mpi_allreduce(mcount0,mcount,1,mpi_integer4,mpi_sum,mpi_comm_world,ierror)
-     pbar4(k)=pbar4(k)/float(mcount)
+     pbar4(k)=pbar4(k)/real(mcount,r_kind)
   end do
 
   psfc015=r015*pbar4(1)
@@ -1160,7 +1160,7 @@ subroutine get_background_glb(mype)
 
   do ilat=1,pf2aP2%nlatf
      do ilon=1,pf2aP2%nlonf
-        if(((float(ilat)-rnf2)**2+(float(ilon)-rnf2)**2)>=rnf212) then
+        if(((real(ilat,r_kind)-rnf2)**2+(real(ilon,r_kind)-rnf2)**2)>=rnf212) then
            p2ilatf(ilat,ilon)=zero
            p3ilatf(ilat,ilon)=zero
         else
@@ -1611,7 +1611,7 @@ subroutine get_aspect_pt(mype)
              cvar=='vp' .or. cvar=='VP' .or. &
              cvar=='t' .or. cvar=='T'
 
-     rk1=float(k1-kthres)
+     rk1=real(k1-kthres,r_kind)
      fblend=half*(one-tanh(rk1))
  
      !--- zonal patch
@@ -1757,7 +1757,7 @@ subroutine get_theta_corrl_lenghts_glb(mype)
      mcount0=lon2*lat2! It's OK to count buffer points
      call mpi_allreduce(pbar4a,pbar4(k),1,mpi_real8,mpi_sum,mpi_comm_world,ierror)
      call mpi_allreduce(mcount0,mcount,1,mpi_integer4,mpi_sum,mpi_comm_world,ierror)
-     pbar4(k)=pbar4(k)/float(mcount)
+     pbar4(k)=pbar4(k)/real(mcount,r_kind)
      call w3fa03(pbar4(k),hgt4(k),tbar4(k),thetabar4(k))
   end do
 
@@ -2605,9 +2605,9 @@ subroutine get_aspect_ens(mype)
 
         nt1=max(1,(nens(k)-1))
  
-        s1=maxval(ensv_p0(:,:,k))/float(nt1)
-        s2=maxval(ensv_p2(:,:,k))/float(nt1)
-        s3=maxval(ensv_p3(:,:,k))/float(nt1)
+        s1=maxval(ensv_p0(:,:,k))/real(nt1,r_kind)
+        s2=maxval(ensv_p2(:,:,k))/real(nt1,r_kind)
+        s3=maxval(ensv_p3(:,:,k))/real(nt1,r_kind)
         smax=max(s1,s2,s3)
 
         if ( nkflag(k)==1 ) then
@@ -3729,13 +3729,13 @@ subroutine ens_intpglb_coeff(iref,jref,enscoeff,mype)
         xg=rlon+one
         yg=rlat+90._r_kind+one
 
-        dxg =xg-float(floor(xg))
-        dyg =yg-float(floor(yg))
+        dxg =xg-real(floor(xg),r_kind)
+        dyg =yg-real(floor(yg),r_kind)
         dxg1=one-dxg
         dyg1=one-dyg
 
-        if (xg >= one .and. xg <= float(jxp) .and.  &
-            yg >= one .and. yg <= float(iy) ) then
+        if (xg >= one .and. xg <= real(jxp,r_kind) .and.  &
+            yg >= one .and. yg <= real(iy,r_kind) ) then
            enscoeff(1,i,j)=dxg1*dyg1
            enscoeff(2,i,j)=dxg1*dyg
            enscoeff(3,i,j)=dxg *dyg1
@@ -3938,7 +3938,7 @@ subroutine ens_uv2psichi(work1,work2)
      vor_s = vor_s + grid_vor( 1,ix)
      vor_n = vor_n + grid_vor(ny,ix)
   end do
-  rnlon = one/float(nlon)
+  rnlon = one/real(nlon,r_kind)
   div_s = div_s*rnlon
   div_n = div_n*rnlon
   vor_s = vor_s*rnlon
diff --git a/src/gsi/atms_spatial_average_mod.f90 b/src/gsi/atms_spatial_average_mod.f90
index b3e4aafc41..639bb8c99c 100644
--- a/src/gsi/atms_spatial_average_mod.f90
+++ b/src/gsi/atms_spatial_average_mod.f90
@@ -841,7 +841,7 @@ SUBROUTINE SFFTCB( X, N, M )
           END DO
           J = J + K
  104  CONTINUE
-      XT = 1.0_r_kind / FLOAT( N )
+      XT = 1.0_r_kind / real( N,r_kind )
       DO 99, I = 1, N
          X(I) = XT * X(I)
  99   CONTINUE
diff --git a/src/gsi/balmod.f90 b/src/gsi/balmod.f90
index 1b9fa9030b..1408530a3f 100644
--- a/src/gsi/balmod.f90
+++ b/src/gsi/balmod.f90
@@ -443,7 +443,7 @@ subroutine prebal_reg(cwcoveqqcov)
              do i=1,lat2
                 l=int(rllat1(i,j))
                 l2=min0(l+1,llmax)
-                dl2=rllat1(i,j)-float(l)
+                dl2=rllat1(i,j)-real(l,r_kind)
                 dl1=one-dl2
                 bvk(i,j,k)=dl1*bvi(l,k)+dl2*bvi(l2,k)
              end do
@@ -465,7 +465,7 @@ subroutine prebal_reg(cwcoveqqcov)
                 do i=1,lat2
                    l=int(rllat1(i,j))
                    l2=min0(l+1,llmax)
-                   dl2=rllat1(i,j)-float(l)
+                   dl2=rllat1(i,j)-real(l,r_kind)
                    dl1=one-dl2
                    agvk(i,j,m,k)=dl1*agvi(l,m,k)+dl2*agvi(l2,m,k)
                 end do
@@ -477,7 +477,7 @@ subroutine prebal_reg(cwcoveqqcov)
              do i=1,lat2
                 l=int(rllat1(i,j))
                 l2=min0(l+1,llmax)
-                dl2=rllat1(i,j)-float(l)
+                dl2=rllat1(i,j)-real(l,r_kind)
                 dl1=one-dl2
                 wgvk(i,j,k)=dl1*wgvi(l,k)+dl2*wgvi(l2,k)
              end do
@@ -972,7 +972,7 @@ subroutine locatelat_reg(mype)
     do j=1,nlon 
        do i=1,nlat   
           if(region_lat(i,j)>=clat_avn(mlat))then
-             rllat(i,j)=float(mlat)
+             rllat(i,j)=real(mlat,r_kind)
              llmax=max0(mlat,llmax)
              llmin=min0(mlat,llmin)
           else if(region_lat(i,j)<clat_avn(1))then
@@ -984,7 +984,7 @@ subroutine locatelat_reg(mype)
                 m1=m+1
                 if((region_lat(i,j)>=clat_avn(m)).and.  &
                    (region_lat(i,j)<clat_avn(m1)))then
-                   rllat(i,j)=float(m)
+                   rllat(i,j)=real(m,r_kind)
                    llmax=max0(m,llmax)
                    llmin=min0(m,llmin)
                    exit
diff --git a/src/gsi/berror.f90 b/src/gsi/berror.f90
index 3a683c9b94..029caffa19 100644
--- a/src/gsi/berror.f90
+++ b/src/gsi/berror.f90
@@ -661,7 +661,7 @@ subroutine init_rftable(mype,rate,nnn,sli,sli1,sli2)
 
     ihwlb=hwlb/tin
     hwlb=ihwlb*tin
-!   tin=(hwle-hwlb)/float(nta-1)
+!   tin=(hwle-hwlb)/real(nta-1,r_kind)
     ntax=(hwle-hwlb)/tin+2
 !   write(6,*)'INIT_RFTABLE:  tin ',ntax,ihwlb,tin,hwlb,hwle
 
@@ -701,7 +701,7 @@ subroutine init_rftable(mype,rate,nnn,sli,sli1,sli2)
        if(iuse(i))then
           nta=nta+1
           ipoint(i)=nta
-          dsh(nta)=one/(float(i-1+ihwlb)*tin)
+          dsh(nta)=one/(real(i-1+ihwlb,r_kind)*tin)
        end if
     end do
 !   write(6,*)'INIT_RFTABLE:  ntax,nta = ',ntax,nta
diff --git a/src/gsi/bkgvar_rewgt.f90 b/src/gsi/bkgvar_rewgt.f90
index ca82882af3..c463723206 100644
--- a/src/gsi/bkgvar_rewgt.f90
+++ b/src/gsi/bkgvar_rewgt.f90
@@ -99,7 +99,7 @@ subroutine bkgvar_rewgt(sfvar,vpvar,tvar,psvar,mype)
   balt   =zero ; bald    =zero ; balps =zero
 
 ! Set count to number of global grid points in quad precision
-  fcount = float(nlat)*float(nlon)
+  fcount = real(nlat*nlon,r_kind)
 
 ! Set parameter for communication
   mm1=mype+1
diff --git a/src/gsi/buddycheck_mod.f90 b/src/gsi/buddycheck_mod.f90
index ae4a29cc63..9c58fd0a44 100644
--- a/src/gsi/buddycheck_mod.f90
+++ b/src/gsi/buddycheck_mod.f90
@@ -200,7 +200,7 @@ subroutine buddy_check_t(is,data,luse,mype,nele,nobs,muse,buddyuse)
   end if
 
 
-  rsig=float(nsig)
+  rsig=real(nsig,r_kind)
   mm1=mype+1
 
   !initialize buddyuse to 1, start by assuming all obs are good!
diff --git a/src/gsi/calc_fov_crosstrk.f90 b/src/gsi/calc_fov_crosstrk.f90
index 75fffb1087..dc9767e850 100644
--- a/src/gsi/calc_fov_crosstrk.f90
+++ b/src/gsi/calc_fov_crosstrk.f90
@@ -812,7 +812,7 @@ subroutine instrument_init(instr, satid, expansion, valid)
  allocate (eccen(1:maxfov(instr)))
  
  do i = 1, npoly
-    psi(i) = two*pi*float(i-1)/float(npoly-1) ! Will connect Npoly points
+    psi(i) = two*pi*real(i-1,r_kind)/real(npoly-1,r_kind) ! Will connect Npoly points
  enddo
  
 ! Precompute angles and sizes for speed. For accurate representation of fov, 
diff --git a/src/gsi/compact_diffs.f90 b/src/gsi/compact_diffs.f90
index ae03407917..0a65e9a515 100644
--- a/src/gsi/compact_diffs.f90
+++ b/src/gsi/compact_diffs.f90
@@ -313,10 +313,10 @@ subroutine stvp2uv(work,idim)
      polsu=polsu+grid3(1 ,ix)*coslon(ix)+grid1(1 ,ix)*sinlon(ix)
      polsv=polsv+grid3(1 ,ix)*sinlon(ix)-grid1(1 ,ix)*coslon(ix)
   end do
-  polnu=polnu/float(nlon)
-  polnv=polnv/float(nlon)
-  polsu=polsu/float(nlon)
-  polsv=polsv/float(nlon)
+  polnu=polnu/real(nlon,r_kind)
+  polnv=polnv/real(nlon,r_kind)
+  polsu=polsu/real(nlon,r_kind)
+  polsv=polsv/real(nlon,r_kind)
 ! work(1 is u, work(2 is v
   do j=1,nlon
      do i=2,nlat-1
@@ -452,7 +452,7 @@ subroutine uv2vordiv(work1,work2)
      vor_s = vor_s + grid_vor( 1,ix)
      vor_n = vor_n + grid_vor(ny,ix)
   end do
-  rnlon = one/float(nlon)
+  rnlon = one/real(nlon,r_kind)
   div_s = div_s*rnlon
   div_n = div_n*rnlon
   vor_s = vor_s*rnlon
@@ -758,10 +758,10 @@ subroutine tstvp2uv(work,idim)
      polnv=polnv+grid3n(ix)*sinlon(ix)+coslon(ix)*grid1n(ix)
      polsv=polsv+grid3s(ix)*sinlon(ix)-coslon(ix)*grid1s(ix)
   end do
-  polnu=polnu/float(nlon)
-  polsu=polsu/float(nlon)
-  polnv=polnv/float(nlon)
-  polsv=polsv/float(nlon)
+  polnu=polnu/real(nlon,r_kind)
+  polsu=polsu/real(nlon,r_kind)
+  polnv=polnv/real(nlon,r_kind)
+  polsv=polsv/real(nlon,r_kind)
   
   do ix=1,nlon
      grid3(ny,ix)=grid3(ny,ix)+polnu*coslon(ix)+polnv*sinlon(ix)
@@ -1357,9 +1357,9 @@ subroutine inisph(r,yor,tau,nx,ny)
 ! Load coefficient array
   ri=one/r
   pih=pi/two
-  pi2onx=pi/float(nxh)
+  pi2onx=pi/real(nxh,r_kind)
   do ix=1,nxh
-     coef(lacoy1+ix-1)=(float(ix)-half)*pi2onx
+     coef(lacoy1+ix-1)=(real(ix,r_kind)-half)*pi2onx
   enddo
 
   call cdcoef(nxh,noq,zero,pi,coef(lacoy1),w&
@@ -1980,10 +1980,10 @@ subroutine compact_dlon(b,dbdx,vector)
         polsu=polsu+grid3(1 ,ix)*coslon(ix)
         polsv=polsv+grid3(1 ,ix)*sinlon(ix)
      end do
-     polnu=two*polnu/float(nlon)
-     polnv=two*polnv/float(nlon)
-     polsu=two*polsu/float(nlon)
-     polsv=two*polsv/float(nlon)
+     polnu=two*polnu/real(nlon,r_kind)
+     polnv=two*polnv/real(nlon,r_kind)
+     polsu=two*polsu/real(nlon,r_kind)
+     polsv=two*polsv/real(nlon,r_kind)
      do ix=1,nlon
         grid3n(ix)= polnu*coslon(ix)+polnv*sinlon(ix)
         grid3s(ix)= polsu*coslon(ix)+polsv*sinlon(ix)
@@ -2073,10 +2073,10 @@ subroutine tcompact_dlon(b,dbdx,vector)
         polsu=polsu+coslon(ix)*grid3s(ix)
         polsv=polsv+sinlon(ix)*grid3s(ix)
      end do
-     polnu=two*polnu/float(nlon)
-     polnv=two*polnv/float(nlon)
-     polsu=two*polsu/float(nlon)
-     polsv=two*polsv/float(nlon)
+     polnu=two*polnu/real(nlon,r_kind)
+     polnv=two*polnv/real(nlon,r_kind)
+     polsu=two*polsu/real(nlon,r_kind)
+     polsv=two*polsv/real(nlon,r_kind)
      do ix=1,nlon
         grid3(ny,ix)=grid3(ny,ix)+coslon(ix)*polnu+sinlon(ix)*polnv
         grid3(1 ,ix)=grid3(1 ,ix)+coslon(ix)*polsu+sinlon(ix)*polsv
@@ -2214,10 +2214,10 @@ subroutine compact_dlat(b,dbdy,vector)
         polsu=polsu+grid4(1 ,ix)*sinlon(ix)
         polsv=polsv-grid4(1 ,ix)*coslon(ix)
      end do
-     polnu=two*polnu/float(nlon)
-     polnv=two*polnv/float(nlon)
-     polsu=two*polsu/float(nlon)
-     polsv=two*polsv/float(nlon)
+     polnu=two*polnu/real(nlon,r_kind)
+     polnv=two*polnv/real(nlon,r_kind)
+     polsu=two*polsu/real(nlon,r_kind)
+     polsv=two*polsv/real(nlon,r_kind)
      do ix=1,nlon
         grid4n(ix)=-polnu*sinlon(ix)+polnv*coslon(ix)
         grid4s(ix)= polsu*sinlon(ix)-polsv*coslon(ix)
@@ -2314,10 +2314,10 @@ subroutine tcompact_dlat(b,dbdy,vector)
         polsu=polsu+sinlon(ix)*grid4s(ix)
         polsv=polsv-coslon(ix)*grid4s(ix)
      end do
-     polnu=two*polnu/float(nlon)
-     polnv=two*polnv/float(nlon)
-     polsu=two*polsu/float(nlon)
-     polsv=two*polsv/float(nlon)
+     polnu=two*polnu/real(nlon,r_kind)
+     polnv=two*polnv/real(nlon,r_kind)
+     polsu=two*polsu/real(nlon,r_kind)
+     polsv=two*polsv/real(nlon,r_kind)
      do ix=1,nlon
         grid4(ny,ix)=grid4(ny,ix)-sinlon(ix)*polnu+coslon(ix)*polnv
         grid4(1 ,ix)=grid4(1 ,ix)+sinlon(ix)*polsu-coslon(ix)*polsv
diff --git a/src/gsi/compute_derived.f90 b/src/gsi/compute_derived.f90
index 46413d9986..32df50445f 100644
--- a/src/gsi/compute_derived.f90
+++ b/src/gsi/compute_derived.f90
@@ -467,7 +467,7 @@ subroutine compute_derived(mype,init_pass)
                     do i=indices%ips,indices%ipe
                        l =max(min(int(rllatf(i,j)),mlat),1)
                        l2=min((l+1),mlat)
-                       dl2=rllatf(i,j)-float(l)
+                       dl2=rllatf(i,j)-real(l,r_kind)
                        dl1=one-dl2
 
                        factk=dl1*corz(l,kvar,nrf3_q)+dl2*corz(l2,kvar,nrf3_q)
diff --git a/src/gsi/compute_qvar3d.f90 b/src/gsi/compute_qvar3d.f90
index 851212d4f5..8babcc82fe 100644
--- a/src/gsi/compute_qvar3d.f90
+++ b/src/gsi/compute_qvar3d.f90
@@ -141,7 +141,7 @@ subroutine compute_qvar3d
               d=20.0_r_kind*rhgues(i,j,k) + one
               n=int(d)
               np=n+1
-              dn2=d-float(n)
+              dn2=d-real(n,r_kind)
               dn1=one-dn2
               n=min0(max(1,n),maxvarq1)
               np=min0(max(1,np),maxvarq1)
@@ -200,7 +200,7 @@ subroutine compute_qvar3d
                  d=-2.0_r_kind*log(cwtmp) + one
                  n=int(d)
                  np=n+1
-                 dn2=d-float(n)
+                 dn2=d-real(n,r_kind)
                  dn1=one-dn2
                  n=min0(max(1,n),30)
                  np=min0(max(1,np),30)
diff --git a/src/gsi/cplr_get_fv3_regional_ensperts.f90 b/src/gsi/cplr_get_fv3_regional_ensperts.f90
index 5b8fb153b6..2645723c62 100644
--- a/src/gsi/cplr_get_fv3_regional_ensperts.f90
+++ b/src/gsi/cplr_get_fv3_regional_ensperts.f90
@@ -646,7 +646,7 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar)
         enddo 
  !
  ! CALCULATE ENSEMBLE MEAN
-        bar_norm = one/float(n_ens_fv3sar)
+        bar_norm = one/real(n_ens_fv3sar,r_kind)
         en_bar(m)%values=en_bar(m)%values*bar_norm
  
  ! Copy pbar to module array.  ps_bar may be needed for vertical localization
@@ -1389,7 +1389,7 @@ subroutine ens_spread_dualres_regional_fv3_regional(this,mype,en_perts,nelen)
             call stop2(9999)
          endif
   
-         sp_norm=(one/float(n_ens))
+         sp_norm=one/real(n_ens,r_kind)
   
          sube%values=zero
   !
diff --git a/src/gsi/cplr_get_pseudo_ensperts.f90 b/src/gsi/cplr_get_pseudo_ensperts.f90
index 9d12165409..f9566c1f45 100644
--- a/src/gsi/cplr_get_pseudo_ensperts.f90
+++ b/src/gsi/cplr_get_pseudo_ensperts.f90
@@ -487,7 +487,7 @@ subroutine get_pseudo_ensperts_wrf(this,en_perts,nelen)
     deallocate(blend)
    
   ! Convert to mean
-    bar_norm = one/float(n_ens)
+    bar_norm = one/real(n_ens,r_kind)
     en_bar%values=en_bar%values*bar_norm
   
     if(write_ens_sprd)then
@@ -1064,15 +1064,15 @@ subroutine get_bgtc_center(bc_lon,bc_lat)
     close(32)
   
     if(sn == 'S')then
-      clat=- one * float(iclat)/10.0_r_kind
+      clat=- one * real(iclat,r_kind)/10.0_r_kind
     else
-      clat=float(iclat)/10.0_r_kind
+      clat=real(iclat,r_kind)/10.0_r_kind
     endif
   
     if(ew == 'W')then
-       clon=360._r_kind - float(iclon)/10.0_r_kind
+       clon=360._r_kind - real(iclon,r_kind)/10.0_r_kind
     else
-       clon=float(iclon)/10.0_r_kind
+       clon=real(iclon,r_kind)/10.0_r_kind
     endif
   
   
diff --git a/src/gsi/cplr_get_wrf_mass_ensperts.f90 b/src/gsi/cplr_get_wrf_mass_ensperts.f90
index 775c9c4195..7b6b3b9aa1 100644
--- a/src/gsi/cplr_get_wrf_mass_ensperts.f90
+++ b/src/gsi/cplr_get_wrf_mass_ensperts.f90
@@ -480,7 +480,7 @@ subroutine get_wrf_mass_ensperts_wrf(this,en_perts,nelen,ps_bar)
 
   !
   ! CALCULATE ENSEMBLE MEAN
-         bar_norm = one/float(n_ens)
+         bar_norm = one/real(n_ens,r_kind)
          en_bar%values=en_bar%values*bar_norm
   
   ! Copy pbar to module array.  ps_bar may be needed for vertical localization
@@ -2189,7 +2189,7 @@ subroutine ens_spread_dualres_regional_wrf(this,mype,en_perts,nelen,en_bar)
             call stop2(999)
          endif
   
-    sp_norm=(one/float(n_ens))
+    sp_norm=one/real(n_ens,r_kind)
   
     sube%values=zero
   !
diff --git a/src/gsi/cplr_get_wrf_nmm_ensperts.f90 b/src/gsi/cplr_get_wrf_nmm_ensperts.f90
index 3e30822639..c9179eedb9 100644
--- a/src/gsi/cplr_get_wrf_nmm_ensperts.f90
+++ b/src/gsi/cplr_get_wrf_nmm_ensperts.f90
@@ -859,7 +859,7 @@ subroutine get_wrf_nmm_ensperts_wrf(this,en_perts,nelen,region_lat_ens,region_lo
   
   !
   ! CALCULATE ENSEMBLE MEAN
-      bar_norm = one/float(n_ens)
+      bar_norm = one/real(n_ens,r_kind)
       en_bar%values=en_bar%values*bar_norm
   
   ! Copy pbar to module array.  ps_bar may be needed for vertical localization
diff --git a/src/gsi/cplr_gfs_ensmod.f90 b/src/gsi/cplr_gfs_ensmod.f90
index 6f7d1184c5..55ac32fce4 100644
--- a/src/gsi/cplr_gfs_ensmod.f90
+++ b/src/gsi/cplr_gfs_ensmod.f90
@@ -661,7 +661,8 @@ subroutine parallel_read_nemsio_state_(en_full,m_cvars2dw,m_cvars3d,nlon,nlat,ns
    enddo
    deallocate(rlons)
 
-   fhour = float(nfhour) + float(nfminute)/r60 + float(nfsecondn)/float(nfsecondd)/r3600
+   fhour = real(nfhour,r_kind) + real(nfminute,r_kind)/r60 + &
+           real(nfsecondn,r_kind)/real(nfsecondd,r_kind)/r3600
    odate(1) = idate(4)  !hour
    odate(2) = idate(2)  !month
    odate(3) = idate(3)  !day
@@ -1107,7 +1108,7 @@ subroutine fillpoles_ss_(temp,nlon,nlat)
       sumn=sumn+temp(nlatm1,i)
       sums=sums+temp(2,i)
    end do
-   rnlon=one/float(nlon)
+   rnlon=one/real(nlon,r_kind)
    sumn_sing=sumn*rnlon
    sums_sing=sums*rnlon
 
@@ -1169,10 +1170,10 @@ subroutine fillpoles_sv_(tempu,tempv,nlon,nlat,clons,slons)
       polsu=polsu+tempu(2,i    )*clons(i)+tempv(2,i    )*slons(i)
       polsv=polsv+tempu(2,i    )*slons(i)-tempv(2,i    )*clons(i)
    end do
-   polnu=polnu/float(nlon)
-   polnv=polnv/float(nlon)
-   polsu=polsu/float(nlon)
-   polsv=polsv/float(nlon)
+   polnu=polnu/real(nlon,r_kind)
+   polnv=polnv/real(nlon,r_kind)
+   polsu=polsu/real(nlon,r_kind)
+   polsv=polsv/real(nlon,r_kind)
    do i=1,nlon
       tempu(nlat,i)= polnu*clons(i)+polnv*slons(i)
       tempv(nlat,i)=-polnu*slons(i)+polnv*clons(i)
diff --git a/src/gsi/cplr_read_wrf_nmm_guess.f90 b/src/gsi/cplr_read_wrf_nmm_guess.f90
index 63b1bac54a..cc1ecf5528 100644
--- a/src/gsi/cplr_read_wrf_nmm_guess.f90
+++ b/src/gsi/cplr_read_wrf_nmm_guess.f90
@@ -2077,7 +2077,7 @@ subroutine read_nems_nmmb_guess_wrf(this,mype)
       !           due to interpolation
        do i=1,lon2
           do j=1,lat2
-             veg_type(j,i,it)=float(nint(veg_type(j,i,it)))
+             veg_type(j,i,it)=real(nint(veg_type(j,i,it)),r_kind)
           end do
        end do
   !                            ! veg frac
@@ -2091,7 +2091,7 @@ subroutine read_nems_nmmb_guess_wrf(this,mype)
       !           due to interpolation
        do i=1,lon2
           do j=1,lat2
-             soil_type(j,i,it)=float(nint(soil_type(j,i,it)))
+             soil_type(j,i,it)=real(nint(soil_type(j,i,it)),r_kind)
           end do
        end do
   
@@ -2102,7 +2102,7 @@ subroutine read_nems_nmmb_guess_wrf(this,mype)
       !           due to interpolation
        do i=1,lon2
           do j=1,lat2
-             smthis(j,i)=float(nint(smthis(j,i)))
+             smthis(j,i)=real(nint(smthis(j,i)),r_kind)
           end do
        end do
   
@@ -2113,7 +2113,7 @@ subroutine read_nems_nmmb_guess_wrf(this,mype)
       !           due to interpolation
        do i=1,lon2
           do j=1,lat2
-             sicethis(j,i)=float(nint(sicethis(j,i)))
+             sicethis(j,i)=real(nint(sicethis(j,i)),r_kind)
           end do
        end do
   
diff --git a/src/gsi/cplr_wrf_netcdf_interface.f90 b/src/gsi/cplr_wrf_netcdf_interface.f90
index 5739b27b73..7c1bc400ec 100644
--- a/src/gsi/cplr_wrf_netcdf_interface.f90
+++ b/src/gsi/cplr_wrf_netcdf_interface.f90
@@ -3388,7 +3388,7 @@ subroutine update_netcdf_mass_wrf(this)
       myname_,'put_att:  START_DATE '//trim(flnm1) )
   call nc_check( nf90_put_att(dh1,nf90_global,'SIMULATION_START_DATE',trim(DateStr1)),&
       myname_,'put_att:  SIMULATION_START_DATE '//trim(flnm1) )
-  call nc_check( nf90_put_att(dh1,nf90_global,'GMT',float(ihour)),&
+  call nc_check( nf90_put_att(dh1,nf90_global,'GMT',real(ihour,r_single)),&
       myname_,'put_att: GMT '//trim(flnm1) )
   call nc_check( nf90_put_att(dh1,nf90_global,'JULYR',iyear),&
       myname_,'put_att: JULYR'//trim(flnm1) )
@@ -3885,7 +3885,7 @@ subroutine update_netcdf_nmm_wrf(this)
       myname_,'put_att:  START_DATE '//trim(adjustl(flnm1)) )
   call nc_check( nf90_put_att(dh1,nf90_global,'SIMULATION_START_DATE',trim(DateStr1)),&
       myname_,'put_att:  SIMULATION_START_DATE '//trim(adjustl(flnm1)) )
-  call nc_check( nf90_put_att(dh1,nf90_global,'GMT',float(ihour)),&
+  call nc_check( nf90_put_att(dh1,nf90_global,'GMT',real(ihour,r_single)),&
       myname_,'put_att: GMT '//trim(adjustl(flnm1)) )
   call nc_check( nf90_put_att(dh1,nf90_global,'JULYR',iyear),&
       myname_,'put_att: JULYR'//trim(adjustl(flnm1)) )
diff --git a/src/gsi/deter_sfc_mod.f90 b/src/gsi/deter_sfc_mod.f90
index e4e77283a4..0e64c9a357 100644
--- a/src/gsi/deter_sfc_mod.f90
+++ b/src/gsi/deter_sfc_mod.f90
@@ -1048,21 +1048,21 @@ subroutine deter_sfc_fov(fov_flag,ifov,instr,ichan,sat_aziang,dlat_earth_deg,&
         exit loop1
      endif
 
-     mid = (float(subgrid_lengths_y)-one)/two + one
-     del = one/ float(subgrid_lengths_y)
+     mid = (real(subgrid_lengths_y,r_kind)-one)/two + one
+     del = one/ real(subgrid_lengths_y,r_kind)
 
      allocate (y_off(subgrid_lengths_y))
 
      do i= 1, subgrid_lengths_y
-        y_off(i) = (float(i)-mid)*del
+        y_off(i) = (real(i,r_kind)-mid)*del
      enddo
 
-     mid = (float(subgrid_lengths_x)-one)/two + one
-     del = one / float(subgrid_lengths_x)
+     mid = (real(subgrid_lengths_x,r_kind)-one)/two + one
+     del = one / real(subgrid_lengths_x,r_kind)
 
      allocate (x_off(subgrid_lengths_x))
      do i= 1, subgrid_lengths_x
-        x_off(i) = (float(i)-mid)*del
+        x_off(i) = (real(i,r_kind)-mid)*del
      enddo
 
 !    Determine the surface characteristics by integrating over the
@@ -1075,9 +1075,9 @@ subroutine deter_sfc_fov(fov_flag,ifov,instr,ichan,sat_aziang,dlat_earth_deg,&
            do i = min_i(j), max_i(j)
               call time_int_sfc(i,j,itsfc,itsfcp,dtsfc,dtsfcp,sfc_mdl)
               do jjj = 1, subgrid_lengths_y
-                 y = float(j) + y_off(jjj)
+                 y = real(j,r_kind) + y_off(jjj)
                  do iii = 1, subgrid_lengths_x
-                    x = float(i) + x_off(iii)
+                    x = real(i,r_kind) + x_off(iii)
                     call txy2ll(x,y,lon_rad,lat_rad)
                     lat_mdl = lat_rad*rad2deg
                     lon_mdl = lon_rad*rad2deg
@@ -1120,7 +1120,7 @@ subroutine deter_sfc_fov(fov_flag,ifov,instr,ichan,sat_aziang,dlat_earth_deg,&
 !              ok here when calculating longitude even if the value is
 !              greater than 360. the ellipse code works from longitude relative
 !              to the center of the fov.
-                    lon_mdl = (float(i)+x_off(iii) - one) * dx_gfs(jj)
+                    lon_mdl = (real(i,r_kind)+x_off(iii) - one) * dx_gfs(jj)
                     if (fov_flag=="crosstrk")then
                        call inside_fov_crosstrk(instr,ifov,sat_aziang, &
                                                dlat_earth_deg,dlon_earth_deg, &
@@ -1916,7 +1916,7 @@ subroutine calc_sfc(sfc_sum,isflg,idomsfc,sfcpct,vfr,sty,vty,sm, &
      vty=zero
   else
      itmp=lbound(sfc_sum%count_vty)-1+maxloc(sfc_sum%count_vty)
-     vty=float(itmp(1))
+     vty=real(itmp(1),r_kind)
   endif
 
 ! soil type is predominate type
@@ -1925,7 +1925,7 @@ subroutine calc_sfc(sfc_sum,isflg,idomsfc,sfcpct,vfr,sty,vty,sm, &
      sty=zero
   else
      itmp=lbound(sfc_sum%count_sty)-1+maxloc(sfc_sum%count_sty)
-     sty=float(itmp(1))
+     sty=real(itmp(1),r_kind)
   endif
 
 ! fields for bare (non-snow covered) land
diff --git a/src/gsi/ens_spread_mod.f90 b/src/gsi/ens_spread_mod.f90
index da70d31ab0..c6528f64ff 100644
--- a/src/gsi/ens_spread_mod.f90
+++ b/src/gsi/ens_spread_mod.f90
@@ -78,7 +78,7 @@ subroutine ens_spread_dualres_regional(en_bar)
           call stop2(999)
        endif
 
-  sp_norm=(one/float(n_ens))
+  sp_norm=(one/real(n_ens,r_kind))
 
   sube%values=zero
 !
diff --git a/src/gsi/general_commvars_mod.f90 b/src/gsi/general_commvars_mod.f90
index daf182fbdd..4304eb6428 100644
--- a/src/gsi/general_commvars_mod.f90
+++ b/src/gsi/general_commvars_mod.f90
@@ -515,7 +515,7 @@ subroutine fill_ns(grid_in,grid_out)
       sumn=sumn+grid_in(i,1)
       sums=sums+grid_in(i,nlatm2)
    end do
-   rnlon=one/float(nlon)
+   rnlon=one/real(nlon,r_kind)
    sumn=sumn*rnlon
    sums=sums*rnlon
 
@@ -620,7 +620,7 @@ subroutine fill2_ns(grid_in,grid_out,nlat,nlon)
       sumn=sumn+grid_in(i,1)
       sums=sums+grid_in(i,nlatm2)
    end do
-   rnlon=one/float(nlon)
+   rnlon=one/real(nlon,r_kind)
    sumn=sumn*rnlon
    sums=sums*rnlon
 
@@ -732,10 +732,10 @@ subroutine filluv_ns(gridu_in,gridv_in,gridu_out,gridv_out)
       polsu=polsu+grid(i,2        )*coslon(i)+grid2(i,2        )*sinlon(i)
       polsv=polsv+grid(i,2        )*sinlon(i)-grid2(i,2        )*coslon(i)
    end do
-   polnu=polnu/float(nlon)
-   polnv=polnv/float(nlon)
-   polsu=polsu/float(nlon)
-   polsv=polsv/float(nlon)
+   polnu=polnu/real(nlon,r_kind)
+   polnv=polnv/real(nlon,r_kind)
+   polsu=polsu/real(nlon,r_kind)
+   polsv=polsv/real(nlon,r_kind)
    do i=1,nlon
       grid (i,nlat)= polnu*coslon(i)+polnv*sinlon(i)
       grid2(i,nlat)=-polnu*sinlon(i)+polnv*coslon(i)
@@ -845,10 +845,10 @@ subroutine filluv2_ns(gridu_in,gridv_in,gridu_out,gridv_out,nlat,nlon,sinlon,cos
       polsu=polsu+grid(i,2        )*coslon(i)+grid2(i,2        )*sinlon(i)
       polsv=polsv+grid(i,2        )*sinlon(i)-grid2(i,2        )*coslon(i)
    end do
-   polnu=polnu/float(nlon)
-   polnv=polnv/float(nlon)
-   polsu=polsu/float(nlon)
-   polsv=polsv/float(nlon)
+   polnu=polnu/real(nlon,r_kind)
+   polnv=polnv/real(nlon,r_kind)
+   polsu=polsu/real(nlon,r_kind)
+   polsv=polsv/real(nlon,r_kind)
    do i=1,nlon
       grid (i,nlat)= polnu*coslon(i)+polnv*sinlon(i)
       grid2(i,nlat)=-polnu*sinlon(i)+polnv*coslon(i)
diff --git a/src/gsi/general_read_fv3atm.f90 b/src/gsi/general_read_fv3atm.f90
index 847d1c4bd3..d3fcfc1b6c 100644
--- a/src/gsi/general_read_fv3atm.f90
+++ b/src/gsi/general_read_fv3atm.f90
@@ -152,7 +152,8 @@ subroutine general_read_fv3atm_nems(grd,sp_a,filename,uvflag,vordivflag,zflag, &
          call stop2(101)
       endif
 
-      fhour = float(nfhour) + float(nfminute)/r60 + float(nfsecondn)/float(nfsecondd)/r3600
+      fhour = real(nfhour,r_kind) + real(nfminute,r_kind)/r60 + &
+              real(nfsecondn,r_kind)/real(nfsecondd,r_kind)/r3600
       odate(1) = idate(4)  !hour
       odate(2) = idate(2)  !month
       odate(3) = idate(3)  !day
diff --git a/src/gsi/general_read_gfsatm.f90 b/src/gsi/general_read_gfsatm.f90
index ffdcd90c79..2971f27148 100755
--- a/src/gsi/general_read_gfsatm.f90
+++ b/src/gsi/general_read_gfsatm.f90
@@ -1227,7 +1227,8 @@ subroutine general_read_gfsatm_nems(grd,sp_a,filename,uvflag,vordivflag,zflag, &
       enddo
       if (mype==0) write(6,*) trim(my_name), ' has_cf = ', has_cf 
 
-      fhour = float(nfhour) + float(nfminute)/r60 + float(nfsecondn)/float(nfsecondd)/r3600
+      fhour = real(nfhour,r_kind) + real(nfminute,r_kind)/r60 + &
+              real(nfsecondn,r_kind)/real(nfsecondd,r_kind)/r3600
       odate(1) = idate(4)  !hour
       odate(2) = idate(2)  !month
       odate(3) = idate(3)  !day
@@ -2094,7 +2095,7 @@ subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag, &
       call read_vardata(filges, 'time', fhour) ! might need to change this to attribute later
                                                ! depends on model changes from
                                                ! Jeff Whitaker
-      fhour = float(nint(fhour))
+      fhour = real(nint(fhour),r_kind)
 
       odate(1) = idate(4)  !hour
       odate(2) = idate(2)  !month
@@ -2918,7 +2919,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
       call read_vardata(filges, 'time', fhour) ! might need to change this to attribute later
                                                ! depends on model changes from
                                                ! Jeff Whitaker
-      fhour = float(nint(fhour))
+      fhour = real(nint(fhour),r_kind)
 
       odate(1) = idate(4)  !hour
       odate(2) = idate(2)  !month
@@ -3813,7 +3814,7 @@ subroutine general_fill_ns(grd,grid_in,grid_out)
       sumn=sumn+grid_in(i,1)
       sums=sums+grid_in(i,nlatm2)
    enddo
-   rnlon=one/float(grd%nlon)
+   rnlon=one/real(grd%nlon,r_kind)
    sumn=sumn*rnlon
    sums=sums*rnlon
 
@@ -3908,10 +3909,10 @@ subroutine general_filluv_ns(grd,slons,clons,gridu_in,gridv_in,gridu_out,gridv_o
       polsu=polsu+gridu_in(i,nlatm2)*clons(i)+gridv_in(i,nlatm2)*slons(i)
       polsv=polsv+gridu_in(i,nlatm2)*slons(i)-gridv_in(i,nlatm2)*clons(i)
    enddo
-   polnu=polnu/float(grd%nlon)
-   polnv=polnv/float(grd%nlon)
-   polsu=polsu/float(grd%nlon)
-   polsv=polsv/float(grd%nlon)
+   polnu=polnu/real(grd%nlon,r_kind)
+   polnv=polnv/real(grd%nlon,r_kind)
+   polsu=polsu/real(grd%nlon,r_kind)
+   polsv=polsv/real(grd%nlon,r_kind)
 
    ! Transfer local work array to output grid
    do k=1,grd%itotsub
@@ -4008,10 +4009,10 @@ subroutine general_fillu_ns(grd,sp,gridu_in,gridv_in,gridu_out)
       polsu=polsu+gridu_in(i,nlatm2)*sp%clons(i)+gridv_in(i,nlatm2)*sp%slons(i)
       polsv=polsv+gridu_in(i,nlatm2)*sp%slons(i)-gridv_in(i,nlatm2)*sp%clons(i)
    enddo
-   polnu=polnu/float(grd%nlon)
-   polnv=polnv/float(grd%nlon)
-   polsu=polsu/float(grd%nlon)
-   polsv=polsv/float(grd%nlon)
+   polnu=polnu/real(grd%nlon,r_kind)
+   polnv=polnv/real(grd%nlon,r_kind)
+   polsu=polsu/real(grd%nlon,r_kind)
+   polsv=polsv/real(grd%nlon,r_kind)
 
    ! Transfer local work array to output grid
    do k=1,grd%itotsub
@@ -4106,10 +4107,10 @@ subroutine general_fillv_ns(grd,sp,gridu_in,gridv_in,gridv_out)
       polsu=polsu+gridu_in(i,nlatm2)*sp%clons(i)+gridv_in(i,nlatm2)*sp%slons(i)
       polsv=polsv+gridu_in(i,nlatm2)*sp%slons(i)-gridv_in(i,nlatm2)*sp%clons(i)
    enddo
-   polnu=polnu/float(grd%nlon)
-   polnv=polnv/float(grd%nlon)
-   polsu=polsu/float(grd%nlon)
-   polsv=polsv/float(grd%nlon)
+   polnu=polnu/real(grd%nlon,r_kind)
+   polnv=polnv/real(grd%nlon,r_kind)
+   polsu=polsu/real(grd%nlon,r_kind)
+   polsv=polsv/real(grd%nlon,r_kind)
 
    ! Transfer local work array to output grid
    do k=1,grd%itotsub
diff --git a/src/gsi/general_read_nemsaero.f90 b/src/gsi/general_read_nemsaero.f90
index 721545bc65..d1e9f525fd 100644
--- a/src/gsi/general_read_nemsaero.f90
+++ b/src/gsi/general_read_nemsaero.f90
@@ -143,7 +143,8 @@ subroutine general_read_nemsaero(grd,sp_a,filename,mype,gfschem_bundle, &
          call stop2(101)
       end if
 
-      fhour = float(nfhour) + float(nfminute)/r60 + float(nfsecondn)/float(nfsecondd)/r3600
+      fhour = real(nfhour,r_kind) + real(nfminute,r_kind)/r60 + & 
+              real(nfsecondn,r_kind)/real(nfsecondd,r_kind)/r3600
       odate(1) = idate(4)  !hour
       odate(2) = idate(2)  !month
       odate(3) = idate(3)  !day
diff --git a/src/gsi/general_spectral_transforms.f90 b/src/gsi/general_spectral_transforms.f90
index 541923e450..76f4c02c42 100644
--- a/src/gsi/general_spectral_transforms.f90
+++ b/src/gsi/general_spectral_transforms.f90
@@ -99,7 +99,7 @@ subroutine general_g2s0_ad(grd,sp,spectral_in,grid_out)
   integer(i_kind) i,j,jj
 
   do i=1,sp%nc
-     spec_work(i)=spectral_in(i)/float(grd%nlon)
+     spec_work(i)=spectral_in(i)/real(grd%nlon,r_kind)
      if(sp%factsml(i))spec_work(i)=zero
   end do
   do i=2*sp%jcap+3,sp%nc
@@ -277,7 +277,7 @@ subroutine general_s2g0_ad(grd,sp,spectral_out,grid_in)
   call general_sptez_s(sp,spec_work,work,-1)
 
   do i=1,sp%nc
-     spec_work(i)=spec_work(i)*float(grd%nlon)
+     spec_work(i)=spec_work(i)*real(grd%nlon,r_kind)
   end do
   do i=2*sp%jcap+3,sp%nc
      spec_work(i)=two*spec_work(i)
@@ -368,7 +368,6 @@ subroutine sfilter(grd,sp,filter,grid)
 
   call general_sptez_s(sp,spec_work,work,-1)
 
-! gnlon=float(grd%nlon)
   gnlon=real(grd%nlon,r_kind)
   do i=1,sp%nc
      spec_work(i)=spec_work(i)*gnlon
@@ -522,8 +521,8 @@ subroutine general_uvg2zds_ad(grd,sp,zsp,dsp,ugrd,vgrd)
   integer(i_kind) i,j,jj
 
   do i=1,sp%nc
-     spcwrk1(i)=dsp(i)/float(grd%nlon)
-     spcwrk2(i)=zsp(i)/float(grd%nlon)
+     spcwrk1(i)=dsp(i)/real(grd%nlon,r_kind)
+     spcwrk2(i)=zsp(i)/real(grd%nlon,r_kind)
      if(sp%factvml(i))then
         spcwrk1(i)=zero
         spcwrk2(i)=zero
@@ -866,8 +865,8 @@ subroutine general_zds2uvg_ad(grd,sp,zsp,dsp,ugrd,vgrd)
   end do
 
   do i=1,sp%nc
-     spcwrk1(i)=spcwrk1(i)*float(grd%nlon)
-     spcwrk2(i)=spcwrk2(i)*float(grd%nlon)
+     spcwrk1(i)=spcwrk1(i)*real(grd%nlon,r_kind)
+     spcwrk2(i)=spcwrk2(i)*real(grd%nlon,r_kind)
   end do
 
   do i=2*sp%jcap+3,sp%nc
diff --git a/src/gsi/general_sub2grid_mod.f90 b/src/gsi/general_sub2grid_mod.f90
index bacea4688b..2285cab005 100644
--- a/src/gsi/general_sub2grid_mod.f90
+++ b/src/gsi/general_sub2grid_mod.f90
@@ -521,7 +521,7 @@ subroutine get_iuse_pe(npe,nz,iuse_pe)
      else                    
         nskip=npe-nz
         if(nskip > 0)then
-          skip2=float(npe)/float(nskip)
+          skip2=real(npe,r_kind)/real(nskip,r_kind)
           point=zero
           do i=1,nskip
             ipoint=min(max(0,nint(point)),npe) 
@@ -878,7 +878,7 @@ subroutine general_deter_subdomain_nolayout(npe,mype,nlat,nlon,regional, &
 !     Compute number of points on full grid and target number of
 !     point per mpi task (pe)
       npts=nlat*nlon
-      anperpe=float(npts)/float(npe)
+      anperpe=real(npts,r_kind)/real(npe,r_kind)
 
 !     Start with square subdomains
       nrnc=sqrt(anperpe)
diff --git a/src/gsi/gengrid_vars.f90 b/src/gsi/gengrid_vars.f90
index a2d352c0b3..adbf510313 100644
--- a/src/gsi/gengrid_vars.f90
+++ b/src/gsi/gengrid_vars.f90
@@ -60,13 +60,13 @@ subroutine gengrid_vars
 ! This is global run, so get global lons, lats, wgtlats, wgtfactlats
 
 ! Set local constants
-     anlon=float(nlon)
+     anlon=real(nlon,r_kind)
      pih=half*pi
      dlon=two*pi/anlon
 
 ! Load grid lat,lon arrays.  rbs2 is used in pcp.
      do i=1,nlon
-        rlons(i)=float(i-1)*dlon
+        rlons(i)=real(i-1,r_kind)*dlon
         coslon(i)=cos(rlons(i))
         sinlon(i)=sin(rlons(i))
      end do
diff --git a/src/gsi/gesinfo.F90 b/src/gsi/gesinfo.F90
index 0aefd34c76..9d287de414 100644
--- a/src/gsi/gesinfo.F90
+++ b/src/gsi/gesinfo.F90
@@ -62,7 +62,8 @@ subroutine gesinfo
 !                                      nfsecondn  FCST Secs (i_kind) numerator
 !                                      nfsecondd  FCST Secs (i_kind) denominator
 !
-!       %fhour = float(nfhour) + float(nfminute)/r60 + float(nfsecondn)/float(nfsecondd)/r3600
+!       %fhour = real(nfhour,r_kind) + real(nfminute,r_kind)/r60 + &
+!                real(nfsecondn,r_kind)/real(nfsecondd,r_kind)/r3600
 !
 ! attributes:
 !   language: f90
@@ -312,8 +313,8 @@ subroutine gesinfo
               nfhour, nfminute, nfsecondn, nfsecondd
            call stop2(99)
         endif
-        gfshead%fhour = float(nfhour) + float(nfminute)/r60 + &
-                        float(nfsecondn)/float(nfsecondd)/r3600
+        gfshead%fhour = real(nfhour,r_kind) + real(nfminute,r_kind)/r60 + &
+                        real(nfsecondn,r_kind)/real(nfsecondd,r_kind)/r3600
 
         gfshead%idate(1) = idate(4)  !hour
         gfshead%idate(3) = idate(3)  !day
@@ -551,7 +552,7 @@ subroutine gesinfo
   ida(:)=0
   jda(:)=0
   fha(:)=zero
-  fha(2)=-float(int(min_offset/60))
+  fha(2)=-real(int(min_offset/60),r_kind)
   fha(3)=-(min_offset+fha(2)*r60)
   ida(1:3)=iadate(1:3)
   ida(5:6)=iadate(4:5)
@@ -582,7 +583,7 @@ subroutine gesinfo
 ! Get time offset
   call time_4dvar(ianldate,time_offset)
 #ifdef RR_CLOUDANALYSIS
-  fha(2)=float(int(min_offset/60))
+  fha(2)=real(int(min_offset/60),r_kind)
   fha(3)=(min_offset-fha(2)*r60)
   time_offset=time_offset+fha(3)/r60
 #endif
diff --git a/src/gsi/get_gefs_ensperts_dualres.f90 b/src/gsi/get_gefs_ensperts_dualres.f90
index bb5ee374af..ca551efa21 100644
--- a/src/gsi/get_gefs_ensperts_dualres.f90
+++ b/src/gsi/get_gefs_ensperts_dualres.f90
@@ -123,7 +123,7 @@ subroutine get_gefs_ensperts_dualres
   im=en_perts(1,1,1)%grid%im
   jm=en_perts(1,1,1)%grid%jm
   km=en_perts(1,1,1)%grid%km
-  bar_norm = one/float(n_ens)
+  bar_norm = one/real(n_ens,r_kind)
   sig_norm=sqrt(one/max(one,n_ens-one))
 
   ! Create temporary communication information for read ensemble routines
@@ -444,7 +444,7 @@ subroutine ens_spread_dualres(en_bar,ibin)
      call stop2(999)
   endif
 
-  sp_norm=(one/float(n_ens))
+  sp_norm=(one/real(n_ens,r_kind))
 
   sube%values=zero
   do n=1,n_ens
diff --git a/src/gsi/get_gefs_for_regional.f90 b/src/gsi/get_gefs_for_regional.f90
index 43a88ef300..cc5e0a2c86 100644
--- a/src/gsi/get_gefs_for_regional.f90
+++ b/src/gsi/get_gefs_for_regional.f90
@@ -304,8 +304,8 @@ subroutine get_gefs_for_regional
      if (nframe /= 0) call error_msg(trim(my_name),trim(filename),'nframe', &
                                      'getfilehead',istop,nframe)
 
-     fhour = float(nfhour) + float(nfminute)/r60 + &
-             float(nfsecondn)/float(nfsecondd)/r3600
+     fhour = real(nfhour,r_kind) + real(nfminute,r_kind)/r60 + &
+             real(nfsecondn,r_kind)/real(nfsecondd,r_kind)/r3600
 
      nlat_gfs=latb+2
      nlon_gfs=lonb
@@ -897,8 +897,8 @@ subroutine get_gefs_for_regional
                  iimin=min(ii,iimin)
                  jjmax=max(jj,jjmax)
                  jjmin=min(jj,jjmin)
-                 dlon_ens=float(jj)
-                 dlat_ens=float(ii)
+                 dlon_ens=real(jj,r_kind)
+                 dlat_ens=real(ii,r_kind)
                  dlon=one+(dlon_ens-one)*ratio_x
                  dlat=one+(dlat_ens-one)*ratio_y
                  call rotate_wind_ll2xy(work_sub(1,i,j,ku),work_sub(1,i,j,kv), &
@@ -992,7 +992,7 @@ subroutine get_gefs_for_regional
   end do
 
 ! Convert to mean
-  bar_norm = one/float(n_ens_gfs)
+  bar_norm = one/real(n_ens_gfs,r_kind)
   do k=1,grd_mix%nsig
      do j=1,grd_mix%lon2
         do i=1,grd_mix%lat2
diff --git a/src/gsi/get_nmmb_ensperts.f90 b/src/gsi/get_nmmb_ensperts.f90
index ece1780c03..4dc3254ccd 100644
--- a/src/gsi/get_nmmb_ensperts.f90
+++ b/src/gsi/get_nmmb_ensperts.f90
@@ -313,7 +313,7 @@ subroutine get_nmmb_ensperts
    end do ! end do over ensemble
 
 ! Convert to mean
-  bar_norm = one/float(n_ens)
+  bar_norm = one/real(n_ens,r_kind)
   en_bar%values=en_bar%values*bar_norm
    
 ! Copy pbar to module array.  ps_bar may be needed for vertical localization
diff --git a/src/gsi/gfs_stratosphere.f90 b/src/gsi/gfs_stratosphere.f90
index 22581b2db0..6014045c76 100644
--- a/src/gsi/gfs_stratosphere.f90
+++ b/src/gsi/gfs_stratosphere.f90
@@ -308,8 +308,8 @@ subroutine mix_gfs_nmmb_vcoords(deta1 ,aeta1 ,eta1 ,deta2 ,aeta2 ,eta2 ,pdtop,pt
       if (nframe /= 0) call error_msg(trim(my_name),trim(filename),'nframe', &
                                          'getfilehead',istop,nframe)
 
-      fhour = float(nfhour) + float(nfminute)/r60 + &
-              float(nfsecondn)/float(nfsecondd)/r3600
+      fhour = real(nfhour,r_kind) + real(nfminute,r_kind)/r60 + &
+              real(nfsecondn,r_kind)/real(nfsecondd,r_kind)/r3600
       write(6,*) ' input filename=',filename
       write(6,*) ' nemsio head: fhour,idate=',fhour,idate
       write(6,*) ' nemsio head: levs=',levs
@@ -1183,8 +1183,8 @@ subroutine add_gfs_stratosphere
          if (  nframe /= 0 ) call error_msg(trim(my_name),trim(filename),'nframe', &
                                             'getfilehead',istop,nframe)
 
-         fhour = float(nfhour) + float(nfminute)/r60 + &
-                 float(nfsecondn)/float(nfsecondd)/r3600
+         fhour = real(nfhour,r_kind) + real(nfminute,r_kind)/r60 + &
+                 real(nfsecondn,r_kind)/real(nfsecondd,r_kind)/r3600
          if ( mype == 0 ) then
             write(6,*) ' input filename=',filename
             write(6,*) ' nemsio head: fhour,idate=',fhour,idate
@@ -1545,8 +1545,8 @@ subroutine add_gfs_stratosphere
             jj=j+grd_mix%jstart(mm1)-2
             ii=min(grd_mix%nlat,max(1,ii))
             jj=min(grd_mix%nlon,max(1,jj))
-            dlon=float(jj)
-            dlat=float(ii)
+            dlon=real(jj,r_kind)
+            dlat=real(ii,r_kind)
             do k=1,nsig_save
                xspli_r(k)=log(prsl_r(i,j,k)*ten)
             enddo
diff --git a/src/gsi/grdcrd.f90 b/src/gsi/grdcrd.f90
index c20e02ce7f..bb655fc68b 100644
--- a/src/gsi/grdcrd.f90
+++ b/src/gsi/grdcrd.f90
@@ -63,7 +63,7 @@ subroutine grdcrd(d,nd,x,nx,flg)
               ix=isrchf(nx-1,x,d(id),flg)-1
            end if
         end if
-        d(id)=float(ix)+(d(id)-x(ix))/(x(ix+1)-x(ix))
+        d(id)=real(ix,r_kind)+(d(id)-x(ix))/(x(ix+1)-x(ix))
      end do
 
 ! Treat special case of nx=1
@@ -135,7 +135,7 @@ subroutine grdcrd1(d,x,nx,flg)
            ix=isrchf(nx-1,x,d,flg)-1
         end if
      end if
-     d=float(ix)+(d-x(ix))/(x(ix+1)-x(ix))
+     d=real(ix,r_kind)+(d-x(ix))/(x(ix+1)-x(ix))
 
 ! Treat special case of nx=1
   elseif (nx==1) then
diff --git a/src/gsi/gsd_update_mod.f90 b/src/gsi/gsd_update_mod.f90
index 7b43f55fb2..35f7663322 100644
--- a/src/gsi/gsd_update_mod.f90
+++ b/src/gsi/gsd_update_mod.f90
@@ -699,7 +699,7 @@ subroutine gsd_gen_coast_prox
                    nip = nip+1
                 end do
                 end do
-                hcoast_prox(1,i,j) = float(nco)/float (nip)
+                hcoast_prox(1,i,j) = real(nco,r_kind)/real(nip,r_kind)
              end if
            end do
         end do
diff --git a/src/gsi/gsdcloudlib_pseudoq_mod.f90 b/src/gsi/gsdcloudlib_pseudoq_mod.f90
index b7544a860c..3b04558da7 100644
--- a/src/gsi/gsdcloudlib_pseudoq_mod.f90
+++ b/src/gsi/gsdcloudlib_pseudoq_mod.f90
@@ -192,9 +192,9 @@ SUBROUTINE cloudCover_Surface_col(mype,nsig, &
                  endif
 
                  ! convert cloud base observation from AGL to ASL
-                 cl_base_ista = float(ocld(6+ic)) + Oelvtn - zh
+                 cl_base_ista = real(ocld(6+ic),r_kind) + Oelvtn - zh
                  if(zh < 1.0_r_kind .and. Oelvtn > 20.0_r_kind &
-                    .and. float(ocld(6+ic)) < 250.0_r_kind) then
+                    .and. real(ocld(6+ic),r_kind) < 250.0_r_kind) then
                     cycle   ! limit the use of METAR station over oceas for low cloud base
                  endif
                  
@@ -267,7 +267,7 @@ SUBROUTINE cloudCover_Surface_col(mype,nsig, &
 ! -- Use visibility for low-level cloud whether
      if (wthr_type < 30   .and. wthr_type > 20 .and. &
          ocld(13)  < 5000 .and. ocld(13) > 1 ) then
-         betav = 3.912_r_kind / (float(ocld(13)) / 1000._r_kind)
+         betav = 3.912_r_kind / (real(ocld(13),r_kind) / 1000._r_kind)
          vis2qc = ( (betav/144.7_r_kind) ** 1.14_r_kind) / 1000._r_kind
      endif  ! cloud or clear
 
diff --git a/src/gsi/guess_grids.F90 b/src/gsi/guess_grids.F90
index bf493a0628..0601959aad 100644
--- a/src/gsi/guess_grids.F90
+++ b/src/gsi/guess_grids.F90
@@ -1783,7 +1783,7 @@ subroutine load_gsdpbl_hgt(mype)
              k=1
              DO while (abs(pbl_height(i,j,jj)) < 0.0001_r_kind)
                if( thetav(k) > thsfc + 1.0_r_kind ) then
-                 pbl_height(i,j,jj) = float(k) - (thetav(k) - (thsfc + 1.0_r_kind))/   &
+                 pbl_height(i,j,jj) = real(k,r_kind) - (thetav(k) - (thsfc + 1.0_r_kind))/   &
                              max((thetav(k)-thetav(k-1)),0.01_r_kind)
                endif
                k=k+1
@@ -2318,7 +2318,7 @@ subroutine guess_grids_stats3d_(name,a,mype)
          end do
       end do
    end do
-   work_a(nsig+1)=float(lon1*lat1)
+   work_a(nsig+1)=real(lon1*lat1,r_kind)
 
    call mpi_allreduce(work_a,work_a1,nsig+1,mpi_rtype,mpi_sum,&
        mpi_comm_world,ierror)
@@ -2386,7 +2386,7 @@ subroutine guess_grids_stats2d_(name,a,mype)
          work_a(1) = work_a(1) + a(i,j)
       end do
    end do
-   work_a(2)=float(lon1*lat1)
+   work_a(2)=real(lon1*lat1,r_kind)
 
    call mpi_allreduce(work_a,work_a1,2,mpi_rtype,mpi_sum,&
        mpi_comm_world,ierror)
diff --git a/src/gsi/hdraobmod.f90 b/src/gsi/hdraobmod.f90
index 00abeb66be..c56b400909 100644
--- a/src/gsi/hdraobmod.f90
+++ b/src/gsi/hdraobmod.f90
@@ -580,7 +580,7 @@ subroutine read_hdraob(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
 !             Add obs reference time, then subtract analysis time to get obs
 !             time relative to analysis
 
-              time_correction=float(minobs-minan)*r60inv
+              time_correction=real(minobs-minan,r_kind)*r60inv
 
            else
               time_correction=zero
diff --git a/src/gsi/hybrid_ensemble_isotropic.F90 b/src/gsi/hybrid_ensemble_isotropic.F90
index fc87026c98..4bad129a72 100644
--- a/src/gsi/hybrid_ensemble_isotropic.F90
+++ b/src/gsi/hybrid_ensemble_isotropic.F90
@@ -247,7 +247,7 @@ subroutine init_rf_z(z_len)
   kap1=rd_over_cp+one
   kapr=one/rd_over_cp
   nxy=grd_ens%latlon11
-  rnsig=float(nsig)
+  rnsig=real(nsig,r_kind)
 
 !    use new factorization:
 
diff --git a/src/gsi/intjcmod.f90 b/src/gsi/intjcmod.f90
index 2b093312ac..f36e9e4b26 100644
--- a/src/gsi/intjcmod.f90
+++ b/src/gsi/intjcmod.f90
@@ -740,7 +740,7 @@ subroutine intjcpdry(rval,sval,nbins,pjc)
   
   it=ntguessig
   mass=zero_quad
-  rcon=(one_quad/(two_quad*float(nlon)))**2
+  rcon=(one_quad/(two_quad*real(nlon,r_quad)))**2
   mm1=mype+1
 
   do n=1,nbins
@@ -1020,7 +1020,7 @@ subroutine intjcpdry2(rval,nbins,mass,pjc)
   integer(i_kind) :: n
   
   it=ntguessig
-  rcon=(one_quad/(two_quad*float(nlon)))**2
+  rcon=(one_quad/(two_quad*real(nlon,r_quad)))**2
   mm1=mype+1
 
   do n=1,nbins
diff --git a/src/gsi/intrp2a.f90 b/src/gsi/intrp2a.f90
index 0d2926b8c7..129b48059f 100644
--- a/src/gsi/intrp2a.f90
+++ b/src/gsi/intrp2a.f90
@@ -56,8 +56,8 @@ subroutine intrp2a(f,g,dx,dy,n,nlevs,mype)
      ix1=int(dx(i))
      iy1=int(dy(i))
      ix1=max(1,min(ix1,nlat))
-     delx=dx(i)-float(ix1)
-     dely=dy(i)-float(iy1)
+     delx=dx(i)-real(ix1,r_kind)
+     dely=dy(i)-real(iy1,r_kind)
      delx=max(zero,min(delx,one))
      ix=ix1-istart(mm1)+2
      iy=iy1-jstart(mm1)+2
@@ -135,8 +135,8 @@ subroutine intrp2a1(f,g,dx,dy,nlevs,mype)
   ix1=int(dx)
   iy1=int(dy)
   ix1=max(1,min(ix1,nlat))
-  delx=dx-float(ix1)
-  dely=dy-float(iy1)
+  delx=dx-real(ix1,r_kind)
+  dely=dy-real(iy1,r_kind)
   delx=max(zero,min(delx,one))
   ix=ix1-istart(mm1)+2
   iy=iy1-jstart(mm1)+2
@@ -211,8 +211,8 @@ subroutine intrp2a11(f,g,dx,dy,mype)
   ix1=int(dx)
   iy1=int(dy)
   ix1=max(1,min(ix1,nlat))
-  delx=dx-float(ix1)
-  dely=dy-float(iy1)
+  delx=dx-real(ix1,r_kind)
+  dely=dy-real(iy1,r_kind)
   delx=max(zero,min(delx,one))
   ix=ix1-istart(mm1)+2
   iy=iy1-jstart(mm1)+2
diff --git a/src/gsi/m_berror_stats.f90 b/src/gsi/m_berror_stats.f90
index 808aee9954..088a7619fe 100644
--- a/src/gsi/m_berror_stats.f90
+++ b/src/gsi/m_berror_stats.f90
@@ -646,7 +646,7 @@ subroutine setcoroz_(coroz,mype)
          enddo
       enddo
    enddo
-   work_oz(nsig+1,mm1)=float(lon1*lat1)
+   work_oz(nsig+1,mm1)=real(lon1*lat1,r_kind)
 
    call mpi_allreduce(work_oz,work_oz1,(nsig+1)*npe,mpi_rtype,mpi_sum,&
       mpi_comm_world,ierror)
@@ -869,7 +869,7 @@ subroutine setcorchem_(cname,corchem,rc)
          enddo
       enddo
    enddo
-   work_chem(nsig+1,mm1)=float(lon1*lat1)
+   work_chem(nsig+1,mm1)=real(lon1*lat1,r_kind)
   
    call mpi_allreduce(work_chem,work_chem1,(nsig+1)*npe,mpi_rtype,mpi_sum,&
         mpi_comm_world,ierror)
diff --git a/src/gsi/m_extOzone.F90 b/src/gsi/m_extOzone.F90
index a28209292f..bf2b137466 100644
--- a/src/gsi/m_extOzone.F90
+++ b/src/gsi/m_extOzone.F90
@@ -1111,7 +1111,7 @@ subroutine ozlev_ncread_(dfile,dtype,ozout,nmrecs,ndata,nodata, gstime,twind)
               ozout(8,ndata)=usage
               ozout(9,ndata)=pob                  ! pressure 
               ozout(10,ndata)=obserr              ! ozone mixing ratio precision in ppmv
-              ozout(11,ndata)=float(ipos(ilev))   ! pointer of obs level index in ozinfo.txt
+              ozout(11,ndata)=real(ipos(ilev),r_kind)   ! pointer of obs level index in ozinfo.txt
               ozout(12,ndata)=levs                ! # of  vertical levels
               ozout(13,ndata)=ppmv                ! ozone mixing ratio in ppmv
            endif
@@ -1424,7 +1424,7 @@ subroutine ozlev_bufrread_(dfile,dtype,dsis, ozout,nmrecs,ndata,nodata, &
            ozout( 8,ndata)=usage1(k)          ! 
            ozout( 9,ndata)=mlspres(k)          ! mls pressure in log(cb)
            ozout(10,ndata)=mlsozpc(k)   ! ozone mixing ratio precision in ppmv
-           ozout(11,ndata)=float(ipos(k))       ! pointer of obs level index in ozinfo.txt
+           ozout(11,ndata)=real(ipos(k),r_kind)       ! pointer of obs level index in ozinfo.txt
            ozout(12,ndata)=nloz         ! # of mls vertical levels
            ozout(13,ndata)=mlsoz(k)     ! ozone mixing ratio in ppmv
         endif
diff --git a/src/gsi/mod_fv3_lola.f90 b/src/gsi/mod_fv3_lola.f90
index 84f8144968..ebe0816c4a 100644
--- a/src/gsi/mod_fv3_lola.f90
+++ b/src/gsi/mod_fv3_lola.f90
@@ -321,10 +321,10 @@ subroutine generate_anl_grid(nx,ny,grid_lon,grid_lont,grid_lat,grid_latt)
 
 !!!!  define analysis A grid  !!!!!!!!!!!!!
   do j=1,nxa
-     xa_a(j)=(float(j-nlonh)-cx)*grid_ratio_fv3_regional
+     xa_a(j)=(real(j-nlonh,r_kind)-cx)*grid_ratio_fv3_regional
   end do
   do i=1,nya
-     ya_a(i)=(float(i-nlath)-cy)*grid_ratio_fv3_regional
+     ya_a(i)=(real(i-nlath,r_kind)-cy)*grid_ratio_fv3_regional
   end do
 
 !!!!!compute fv3 to A grid interpolation parameters !!!!!!!!!
diff --git a/src/gsi/mod_vtrans.f90 b/src/gsi/mod_vtrans.f90
index 3c7d8af8f4..f9538e0735 100644
--- a/src/gsi/mod_vtrans.f90
+++ b/src/gsi/mod_vtrans.f90
@@ -269,7 +269,7 @@ subroutine create_vtrans(mype)
 
 !   count:
 !  Not clear if area weighting would be better.
-    count=one/float(nlat*nlon)
+    count=one/real(nlat*nlon,r_kind)
 
     ier=0
     call gsi_bundlegetpointer (gsi_metguess_bundle(ntguessig),'ps',ges_ps_nt, istatus)
diff --git a/src/gsi/mod_wrfmass_to_a.f90 b/src/gsi/mod_wrfmass_to_a.f90
index 716ff37f77..594f128709 100644
--- a/src/gsi/mod_wrfmass_to_a.f90
+++ b/src/gsi/mod_wrfmass_to_a.f90
@@ -312,7 +312,7 @@ subroutine wrfmass_obs_to_a8(obsba,nreal,maxobs,ilat,ilon,numobs)
      i=int(ria(n)+0.5_r_kind)
      j=int(rja(n)+0.5_r_kind)
 
-     adist=(ria(n)-float(i))*(ria(n)-float(i))+(rja(n)-float(j))*(rja(n)-float(j))
+     adist=(ria(n)-real(i,r_kind))**2+(rja(n)-real(j,r_kind))**2
      if(adist < dist(i,j)) then
         dist(i,j)=adist
         ija(i,j)=n
@@ -324,8 +324,8 @@ subroutine wrfmass_obs_to_a8(obsba,nreal,maxobs,ilat,ilon,numobs)
      do i=1,nxa
         if(ija(i,j) > 0) then
            n=n+1
-           obsba(ilon,n)=float(i)
-           obsba(ilat,n)=float(j)
+           obsba(ilon,n)=real(i,r_kind)
+           obsba(ilat,n)=real(j,r_kind)
            do k=3,nreal
               obsba(k,n)=obsa(k,ija(i,j))
            enddo
diff --git a/src/gsi/mp_compact_diffs_mod1.f90 b/src/gsi/mp_compact_diffs_mod1.f90
index bb924441b7..292986a56a 100644
--- a/src/gsi/mp_compact_diffs_mod1.f90
+++ b/src/gsi/mp_compact_diffs_mod1.f90
@@ -1550,8 +1550,8 @@ subroutine mp_compact_dlon(b,dbdx,vector)
                     polu=polu+grid3(ix)*coslon(ix)
                     polv=polv+grid3(ix)*sinlon(ix)
                  end do
-                 polu=polu/float(nlon)
-                 polv=polv/float(nlon)
+                 polu=polu/real(nlon,r_kind)
+                 polv=polv/real(nlon,r_kind)
                  do ix=1,nlon
                     grid3pol(ix)=polu*coslon(ix)+polv*sinlon(ix)
                  end do
@@ -1673,8 +1673,8 @@ subroutine mp_compact_dlon_ad(b,dbdx,vector)
                     polu=polu+grid3pol(ix)*coslon(ix)
                     polv=polv+grid3pol(ix)*sinlon(ix)
                  end do
-                 polu=polu/float(nlon)
-                 polv=polv/float(nlon)
+                 polu=polu/real(nlon,r_kind)
+                 polv=polv/real(nlon,r_kind)
                  do ix=1,nlon
                     grid3(ix)=grid3(ix)+polu*coslon(ix)+polv*sinlon(ix)
                  end do
@@ -1977,8 +1977,8 @@ subroutine mp_uv_pole(u,v)
            polsu=polsu+u(2,ix,k)*coslon(ix)+v(2,ix,k)*sinlon(ix)
            polsv=polsv+u(2,ix,k)*sinlon(ix)-v(2,ix,k)*coslon(ix)
         end do
-        polsu=polsu/float(nlon)
-        polsv=polsv/float(nlon)
+        polsu=polsu/real(nlon,r_kind)
+        polsv=polsv/real(nlon,r_kind)
         do ix=1,nlon
            u(1,ix,k)=polsu*coslon(ix)+polsv*sinlon(ix)
            v(1,ix,k)=polsu*sinlon(ix)-polsv*coslon(ix)
@@ -1993,8 +1993,8 @@ subroutine mp_uv_pole(u,v)
            polnu=polnu+u(1,ix,k)*coslon(ix)-v(1,ix,k)*sinlon(ix)
            polnv=polnv+u(1,ix,k)*sinlon(ix)+v(1,ix,k)*coslon(ix)
         end do
-        polnu=polnu/float(nlon)
-        polnv=polnv/float(nlon)
+        polnu=polnu/real(nlon,r_kind)
+        polnv=polnv/real(nlon,r_kind)
         do ix=1,nlon
            u(2,ix,k)= polnu*coslon(ix)+polnv*sinlon(ix)
            v(2,ix,k)=-polnu*sinlon(ix)+polnv*coslon(ix)
@@ -2055,8 +2055,8 @@ subroutine mp_uv_pole_ad(u,v)
            u(1,ix,k)=zero
            v(1,ix,k)=zero
         end do
-        polsu=polsu/float(nlon)
-        polsv=polsv/float(nlon)
+        polsu=polsu/real(nlon,r_kind)
+        polsv=polsv/real(nlon,r_kind)
         do ix=1,nlon
            u(2,ix,k)=u(2,ix,k)+polsu*coslon(ix)+polsv*sinlon(ix)
            v(2,ix,k)=v(2,ix,k)+polsu*sinlon(ix)-polsv*coslon(ix)
@@ -2073,8 +2073,8 @@ subroutine mp_uv_pole_ad(u,v)
            u(2,ix,k)=zero
            v(2,ix,k)=zero
         end do
-        polnu=polnu/float(nlon)
-        polnv=polnv/float(nlon)
+        polnu=polnu/real(nlon,r_kind)
+        polnv=polnv/real(nlon,r_kind)
         do ix=1,nlon
            u(1,ix,k)=u(1,ix,k)+polnu*coslon(ix)+polnv*sinlon(ix)
            v(1,ix,k)=v(1,ix,k)-polnu*sinlon(ix)+polnv*coslon(ix)
diff --git a/src/gsi/ncepgfs_ghg.f90 b/src/gsi/ncepgfs_ghg.f90
index 6c5fa7bb9d..9b34d8aa6e 100644
--- a/src/gsi/ncepgfs_ghg.f90
+++ b/src/gsi/ncepgfs_ghg.f90
@@ -326,7 +326,7 @@ subroutine read_gfsco2  &
             do i=1,nmxlon
                co2diff= co2_sav2(i,j,k)-co2_sav1(i,j,k)
                co2rate= co2diff/ndmax
-               co2_Tintrp(i,j,k)= co2_sav1(i,j,k)+ co2rate*float(idd-1)
+               co2_Tintrp(i,j,k)= co2_sav1(i,j,k)+ co2rate*real(idd-1,r_kind)
             enddo
          enddo
       enddo
@@ -558,7 +558,7 @@ subroutine read_ch4n2oco &
             do i=1,nmaxlon
                ghgdiff= ghg_sav2(1,j,k)-ghg_sav1(1,j,k)
                ghgrate= ghgdiff/ndmax
-               ghg_Tintrp(1,j,k)= ghg_sav1(1,j,k)+ ghgrate*float(idd-1)
+               ghg_Tintrp(1,j,k)= ghg_sav1(1,j,k)+ ghgrate*real(idd-1,r_kind)
             enddo
          enddo
       enddo
diff --git a/src/gsi/ncepgfs_io.f90 b/src/gsi/ncepgfs_io.f90
index 59be6d3925..dd46916039 100644
--- a/src/gsi/ncepgfs_io.f90
+++ b/src/gsi/ncepgfs_io.f90
@@ -2640,9 +2640,9 @@ subroutine write_ens_sfc_nst(mype_so,dsfct)
                 allocate(slatx(jmax),wlatx(jmax))
                 allocate(rlats_ens_sfc(nlat_ens_sfc),rlons_ens_sfc(nlon_ens_sfc))
                 call splat(4,jmax,slatx,wlatx)
-                dlon=two*pi/float(nlon_ens_sfc)
+                dlon=two*pi/real(nlon_ens_sfc,r_kind)
                 do i=1,nlon_ens_sfc
-                   rlons_ens_sfc(i)=float(i-1)*dlon
+                   rlons_ens_sfc(i)=real(i-1,r_kind)*dlon
                 end do
                 do i=1,(nlat_ens_sfc-1)/2
                    rlats_ens_sfc(i+1)=-asin(slatx(i))
@@ -3074,9 +3074,9 @@ subroutine write_ens_dsfct(mype_so,dsfct)
           allocate(slatx(jmax),wlatx(jmax))
           allocate(rlats_ens_sfc(nlat_ens_sfc),rlons_ens_sfc(nlon_ens_sfc))
           call splat(4,jmax,slatx,wlatx)
-          dlon=two*pi/float(nlon_ens_sfc)
+          dlon=two*pi/real(nlon_ens_sfc,r_kind)
           do i=1,nlon_ens_sfc
-             rlons_ens_sfc(i)=float(i-1)*dlon
+             rlons_ens_sfc(i)=real(i-1,r_kind)*dlon
           end do
           do i=1,(nlat_ens_sfc-1)/2
              rlats_ens_sfc(i+1)=-asin(slatx(i))
@@ -3247,7 +3247,7 @@ subroutine glbave(fld,ave)
         enddo
      enddo
   enddo
-  xave=xave/(two_quad*float(nlon))
+  xave=xave/(two_quad*real(nlon,r_quad))
   call mpl_allreduce(size(ave,1),qpvals=xave)
   ave=xave
   deallocate(xave)
diff --git a/src/gsi/ncepnems_io.f90 b/src/gsi/ncepnems_io.f90
index 73f06ffd07..595f07e152 100755
--- a/src/gsi/ncepnems_io.f90
+++ b/src/gsi/ncepnems_io.f90
@@ -82,7 +82,8 @@ module ncepnems_io
 !                                      nfsecondn  FCST Secs (i_kind) numerator
 !                                      nfsecondd  FCST Secs (i_kind) denominator
 !
-!       %fhour = float(nfhour) + float(nfminute)/r60 + float(nfsecondn)/float(nfsecondd)/r3600
+!       %fhour = real(nfhour,r_kind) + real(nfminute,r_kind)/r60 + &
+!                real(nfsecondn,r_kind)/real(nfsecondd,r_kind)/r3600
 !
 !   nframe     - nframe is the number of grids extend outward from the
 !                edge of modeling domain.
@@ -835,7 +836,8 @@ subroutine read_atm_ (grd,filename,sp_a,uvflag,vordivflag,zflag, &
          call stop2(101)
       end if
 
-      fhour = float(nfhour) + float(nfminute)/r60 + float(nfsecondn)/float(nfsecondd)/r3600
+      fhour = real(nfhour,r_kind) + real(nfminute,r_kind)/r60 + &
+            real(nfsecondn,r_kind)/real(nfsecondd,r_kind)/r3600
       odate(1) = idate(4)  !hour
       odate(2) = idate(2)  !month
       odate(3) = idate(3)  !day
@@ -1300,7 +1302,8 @@ subroutine read_sfc_(sfct,soil_moi,sno,soil_temp,veg_frac,fact10,sfc_rough, &
           call stop2(102)
        end if
 
-       fhour = float(nfhour) + float(nfminute)/r60 + float(nfsecondn)/float(nfsecondd)/r3600
+       fhour = real(nfhour,r_kind) + real(nfminute,r_kind)/r60 + &
+               real(nfsecondn,r_kind)/real(nfsecondd,r_kind)/r3600
        odate(1) = idate(4)  !hour
        odate(2) = idate(2)  !month
        odate(3) = idate(3)  !day
@@ -1699,7 +1702,8 @@ subroutine read_sfc_anl_(isli_anl)
        call stop2(102)
     end if
 
-    fhour = float(nfhour) + float(nfminute)/r60 + float(nfsecondn)/float(nfsecondd)/r3600
+    fhour = real(nfhour,r_kind) + real(nfminute,r_kind)/r60 + &
+            real(nfsecondn,r_kind)/real(nfsecondd,r_kind)/r3600
     odate(1) = idate(4)  !hour
     odate(2) = idate(2)  !month
     odate(3) = idate(3)  !day
@@ -2029,7 +2033,8 @@ subroutine read_nst_ (tref,dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d)
           call stop2(istop)
        end if
 
-       fhour = float(nfhour) + float(nfminute)/r60 + float(nfsecondn)/float(nfsecondd)/r3600
+       fhour = real(nfhour,r_kind) + real(nfminute,r_kind)/r60 + &
+               real(nfsecondn,r_kind)/real(nfsecondd,r_kind)/r3600
        odate(1) = idate(4)  !hour
        odate(2) = idate(2)  !month
        odate(3) = idate(3)  !day
@@ -5469,8 +5474,8 @@ subroutine tran_gfssfc(ain,aout,lonb,latb)
      sumn = ain(i,1)    + sumn
      sums = ain(i,latb) + sums
   end do
-  sumn = sumn/float(lonb)
-  sums = sums/float(lonb)
+  sumn = sumn/real(lonb,r_kind)
+  sums = sums/real(lonb,r_kind)
 !  Transfer from local work array to surface guess array
   do j = 1,lonb
      aout(1,j)=sums
diff --git a/src/gsi/netcdfgfs_io.f90 b/src/gsi/netcdfgfs_io.f90
index 41e8f33e03..8aa6c4cc6a 100644
--- a/src/gsi/netcdfgfs_io.f90
+++ b/src/gsi/netcdfgfs_io.f90
@@ -3202,8 +3202,8 @@ subroutine tran_gfsncsfc(ain,aout,lonb,latb)
      sumn = ain(i,1)    + sumn
      sums = ain(i,latb) + sums
   end do
-  sumn = sumn/float(lonb)
-  sums = sums/float(lonb)
+  sumn = sumn/real(lonb,r_kind)
+  sums = sums/real(lonb,r_kind)
 !  Transfer from local work array to surface guess array
   do j = 1,lonb
      aout(1,j)=sums
diff --git a/src/gsi/nlmsas_ad.f90 b/src/gsi/nlmsas_ad.f90
index be8213b340..7ad413aa6b 100644
--- a/src/gsi/nlmsas_ad.f90
+++ b/src/gsi/nlmsas_ad.f90
@@ -707,10 +707,10 @@ subroutine nlmsas_ad_im_ix_(im,ix,km,jcap,delt,del,sl,rcs,&
   pdpdwn  = zero
   pdetrn  = 200._r_kind
   xlambu  = 1.e-4_r_kind
-  fjcap = (float(jcap) / 126._r_kind) ** 2
+  fjcap = (real(jcap,r_kind) / 126._r_kind) ** 2
   val   =  one
   fjcap = max(fjcap,val)
-  fkm = (float(km) / 28._r_kind) ** 2
+  fkm = (real(km,r_kind) / 28._r_kind) ** 2
   fkm = max(fkm,one)
   w1l = -8.e-3_r_kind
   w2l = -4.e-2_r_kind
@@ -1058,7 +1058,7 @@ subroutine nlmsas_ad_im_ix_(im,ix,km,jcap,delt,del,sl,rcs,&
 
 ! Select cloud from ensemble
   do i=1,im
-     kt2(i) = nint(xkt2(i)*float(ktcon(i)-jmin(i))-half) + jmin(i) + 1
+     kt2(i) = nint(xkt2(i)*real(ktcon(i)-jmin(i),r_kind)-half) + jmin(i) + 1
      tem1 = hcko(jmin(i),i) - hesol(kt2(i),i)
      tem2 = sumz(kt2(i),i) * hesol(kt2(i),i) - sumh(kt2(i),i)
      if (abs(tem2)  >  0.000001_r_kind) then
diff --git a/src/gsi/polcarf.f90 b/src/gsi/polcarf.f90
index 4ed450daa1..f2038d9d95 100644
--- a/src/gsi/polcarf.f90
+++ b/src/gsi/polcarf.f90
@@ -617,7 +617,7 @@ subroutine polcas(afg,axr,nxem,norm,naxr,wtaxs,wtxrs,inaxs,inxrs,nf,mr,nr)
   do i=0,naxr-1
      valp=valp+axr(i,mr+1)
   end do
-  valp=valp/float(naxr)
+  valp=valp/real(naxr,r_kind)
   do i=0,naxr-1
      axr(i,mr)=valp
   end do
@@ -692,7 +692,7 @@ subroutine polcasa(afg,axr,nxem,norm,naxr,wtaxs,wtxrs,inaxs,inxrs,nf,mr,nr)
   do i=0,naxr-1
      valp=valp+axr(i,mr)
   end do
-  valp=valp/float(naxr)
+  valp=valp/real(naxr,r_kind)
   do i=0,naxr-1
      axr(i,mr)=zero
      axr(i,mr+1)=axr(i,mr+1)+valp
diff --git a/src/gsi/prewgt.f90 b/src/gsi/prewgt.f90
index 38d61050fe..7f3d2fc091 100644
--- a/src/gsi/prewgt.f90
+++ b/src/gsi/prewgt.f90
@@ -273,7 +273,7 @@ subroutine prewgt(mype)
      end do
      do j=1,lon2
         do i=1,lat2
-           temp(i,j)=float(isli2(i,j))
+           temp(i,j)=real(isli2(i,j),r_kind)
         end do
      end do
 
@@ -585,7 +585,7 @@ subroutine prewgt(mype)
 ! rearth_equator is the equatorial radius from a 1999 IAG report.  The
 ! horizontal scales are defined at the equator, hence the need for the
 ! equatorial radius.
-  s2u=(two*pi*rearth_equator)/float(nlon)
+  s2u=(two*pi*rearth_equator)/real(nlon,r_kind)
 
 
   allocate(sli(ny,nx,2,nnnn1o),sli1(-nf:nf,-nf:nf,2,nnnn1o), &
diff --git a/src/gsi/prewgt_reg.f90 b/src/gsi/prewgt_reg.f90
index bde1c9dece..1da89a703b 100644
--- a/src/gsi/prewgt_reg.f90
+++ b/src/gsi/prewgt_reg.f90
@@ -452,7 +452,7 @@ subroutine prewgt_reg(mype)
               d=region_lat(il,jl)*rad2deg+90._r_kind
               l=int(d)
               l2=l+1
-              dl2=d-float(l)
+              dl2=d-real(l,r_kind)
               dl1=one-dl2
               do k=1,nsig
                  dssv(i,j,k,n)=(dl1*ozmzt(l,k)+dl2*ozmzt(l2,k))*dsv(1,k,llmin)
@@ -581,7 +581,7 @@ subroutine prewgt_reg(mype)
                 do i=1,lon2
                     l=int(rllat1(j,i))
                     l2=min0(l+1,llmax)
-                    dl2=rllat1(j,i)-float(l)
+                    dl2=rllat1(j,i)-real(l,r_kind)
                     dl1=one-dl2
                     do k=1,nsig
                         dssv(j,i,k,n)=dl1*dsv(i,k,l)+dl2*dsv(i,k,l2)
@@ -604,7 +604,7 @@ subroutine prewgt_reg(mype)
            do i=1,lon2
               l=int(rllat1(j,i))
               l2=min0(l+1,llmax)
-              dl2=rllat1(j,i)-float(l)
+              dl2=rllat1(j,i)-real(l,r_kind)
               dl1=one-dl2
               do k=1,nsig
                  dssv(j,i,k,n)=dl1*dsv(i,k,l)+dl2*dsv(i,k,l2)
@@ -662,7 +662,7 @@ subroutine prewgt_reg(mype)
         do i=1,lon2
            l=int(rllat1(j,i))
            l2=min0(l+1,llmax)
-           dl2=rllat1(j,i)-float(l)
+           dl2=rllat1(j,i)-real(l,r_kind)
            dl1=one-dl2
            dssvs(j,i,n)=dl1*dsvs(i,l)+dl2*dsvs(i,l2)
            if (mvars>=2.and.n==nrf2_sst) then
@@ -738,7 +738,7 @@ subroutine prewgt_reg(mype)
                     do j=1,ny
                        l=int(rllat(j,i))
                        lp=min0(l+1,llmax)
-                       dl2=rllat(j,i)-float(l)
+                       dl2=rllat(j,i)-real(l,r_kind)
                        dl1=one-dl2
                        fact=one/(dl1*hwll(l,k1,nn)+dl2*hwll(lp,k1,nn))
                        slw((i-1)*ny+j,k)=slw((i-1)*ny+j,1)*fact**2
@@ -778,7 +778,7 @@ subroutine prewgt_reg(mype)
                  do j=1,ny
                     l=int(rllat(j,i))
                     lp=min0(l+1,llmax)
-                    dl2=rllat(j,i)-float(l)
+                    dl2=rllat(j,i)-real(l,r_kind)
                     dl1=one-dl2
                     fact=cc/(dl1*hwllp(l,nn)+dl2*hwllp(lp,nn))
                     slw((i-1)*ny+j,k)=slw((i-1)*ny+j,1)*fact**2
diff --git a/src/gsi/q_diag.f90 b/src/gsi/q_diag.f90
index 15ef49c6b5..fe5f2f0e94 100644
--- a/src/gsi/q_diag.f90
+++ b/src/gsi/q_diag.f90
@@ -147,7 +147,7 @@ subroutine q_diag(it,mype)
      call load_grid(work_pw,grid_pw)
      globps=zero
      globpw=zero
-     rlon=one/float(nlon)
+     rlon=one/real(nlon,r_kind)
      do jj=2,nlat-1
         j=jj-1
         fmeanps=zero
diff --git a/src/gsi/radinfo.f90 b/src/gsi/radinfo.f90
index ffc4641696..ede58b9bca 100644
--- a/src/gsi/radinfo.f90
+++ b/src/gsi/radinfo.f90
@@ -1454,7 +1454,7 @@ real(r_kind) function rnad_pos(isis,iscan,jch)
 
       piece=-0.625_r_kind
       if (mod(iscan,2) == 1) piece = 0.625_r_kind
-      rnad_pos=radstart(jch)+radstep(jch)*float((iscan-1)/2)+piece
+      rnad_pos=radstart(jch)+radstep(jch)*real((iscan-1)/2,r_kind)+piece
 
    else
 
@@ -1466,7 +1466,7 @@ real(r_kind) function rnad_pos(isis,iscan,jch)
       else
          ifov=iscan
       end if
-      rnad_pos=radstart(jch)+radstep(jch)*float(ifov-1)
+      rnad_pos=radstart(jch)+radstep(jch)*real(ifov-1,r_kind)
 
    end if
 
@@ -2034,7 +2034,7 @@ subroutine init_predx
                         tlap2(jj) = tlap0(jj) + tlap1(jj)/tsum(jj)
                         count_tlapmean(jj)=count_tlapmean(jj)+one
                      elseif (tcnt(jj)>0) then
-                        ratio = max(zero,min(tcnt(jj)/float(nthreshold),one))
+                        ratio = max(zero,min(tcnt(jj)/real(nthreshold,r_kind),one))
                         tsum(jj)=ratio*tsum(jj)+tsum0(jj)
 !                       tlap2(jj) = tlap0(jj) + ratio*wgtlap*tlap1(jj)/tsum(jj)
                         tlap2(jj) = tlap0(jj) + ratio*tlap1(jj)/tsum(jj)
diff --git a/src/gsi/raflib.f90 b/src/gsi/raflib.f90
index 227bda3cb1..eb6d399d21 100644
--- a/src/gsi/raflib.f90
+++ b/src/gsi/raflib.f90
@@ -4488,14 +4488,14 @@ SUBROUTINE EIGEN(A,R,N,MV)
       end do
       if(anorm>zero) then
          ANORM=1.414_r_kind*SQRT(ANORM)
-         ANRMX=ANORM*RANGE/FLOAT(N)
+         ANRMX=ANORM*RANGE/real(N,r_kind)
 !
 !        INITIALIZE INDICATORS AND COMPUTE THRESHOLD, THR
 !
          IND=0
          THR=ANORM
          loop1: do
-            THR=THR/FLOAT(N)
+            THR=THR/real(N,r_kind)
             loop2: do
                L=1
                loop3: do
diff --git a/src/gsi/rdgrbsst.f90 b/src/gsi/rdgrbsst.f90
index 8be85ad608..29e5346eae 100644
--- a/src/gsi/rdgrbsst.f90
+++ b/src/gsi/rdgrbsst.f90
@@ -132,14 +132,14 @@ subroutine rdgrbsst(file_sst,mlat_sst,mlon_sst,&
 
 ! Get lat_sst & lon_sst
   do i = 2, nlat_sst - 1
-     rlats_sst(i) = (xsst0 + float(i-2)*dres)*deg2rad
+     rlats_sst(i) = (xsst0 + real(i-2,r_kind)*dres)*deg2rad
   enddo
 
   rlats_sst(1)        = -90.0_r_kind*deg2rad
   rlats_sst(nlat_sst) =  90.0_r_kind*deg2rad
   
   do j = 2, nlon_sst - 1
-     rlons_sst(j) = (ysst0 + float(j-2)*dres)*deg2rad
+     rlons_sst(j) = (ysst0 + real(j-2,r_kind)*dres)*deg2rad
   enddo
 
   rlons_sst(1)        = -half*dres*deg2rad                              ! 1
@@ -184,8 +184,8 @@ subroutine rdgrbsst(file_sst,mlat_sst,mlon_sst,&
      sums = sums + sst(j,2)
      sumn = sumn + sst(j,nlat_sst-1)
   end do
-  sums = sums / float(i)
-  sumn = sumn / float(i)
+  sums = sums / real(i,r_kind)
+  sumn = sumn / real(i,r_kind)
   do j = 2,nlon_sst-1
      sst(j,1)        = sums
      sst(j,nlat_sst) = sumn
diff --git a/src/gsi/read_airs.f90 b/src/gsi/read_airs.f90
index 16e890abd1..8dd22ffd5e 100644
--- a/src/gsi/read_airs.f90
+++ b/src/gsi/read_airs.f90
@@ -790,7 +790,7 @@ subroutine read_airs(mype,val_airs,ithin,isfcalc,rmesh,jsatid,gstime,&
         endif
 
         sol_aziang = aquaspot(2)
-        lza = (start + float(ifov-1)*step)*deg2rad
+        lza = (start + real(ifov-1,r_kind)*step)*deg2rad
 !
 !       interpolate NSST variables to Obs. location and get dtw, dtc, tz_tr
 !
diff --git a/src/gsi/read_amsr2.f90 b/src/gsi/read_amsr2.f90
index a7b27abccc..9d8d4944d9 100644
--- a/src/gsi/read_amsr2.f90
+++ b/src/gsi/read_amsr2.f90
@@ -566,7 +566,7 @@ subroutine read_amsr2(mype,val_amsr2,ithin,rmesh,jsatid,gstime,&
        if(.not. regional .and. dist1 > 0.75_r_kind) cycle obsloop
     endif
 
-    crit1 = crit1 + 10._r_kind * float(iskip)
+    crit1 = crit1 + 10._r_kind * real(iskip,r_kind)
     call checkob(dist1,crit1,itx,iuse)
     if(.not. iuse) then
        cycle obsloop
diff --git a/src/gsi/read_atms.f90 b/src/gsi/read_atms.f90
index 47b675b3a3..c6ed159068 100644
--- a/src/gsi/read_atms.f90
+++ b/src/gsi/read_atms.f90
@@ -401,7 +401,7 @@ subroutine read_atms(mype,val_tovs,ithin,isfcalc,&
 
 !          inflate selection value for ears_db data
            crit0 = 0.01_r_kind
-           if ( llll > 1 ) crit0 = crit0 + r100 * float(llll)
+           if ( llll > 1 ) crit0 = crit0 + r100 * real(llll,r_kind)
 
            call ufbint(lnbufr,bfr1bhdr,n1bhdr,1,iret,hdr1b)
 
@@ -455,7 +455,7 @@ subroutine read_atms(mype,val_tovs,ithin,isfcalc,&
            lza = bfr2bhdr(1)*deg2rad      ! local zenith angle
            if(ifov <= 48)    lza=-lza
 
-           panglr=(start+float(ifov-1)*step)*deg2rad
+           panglr=(start+real(ifov-1,r_kind)*step)*deg2rad
            satellite_height=bfr1bhdr(13)
 !          Ensure orbit height is reasonable
            if (satellite_height < 780000.0_r_kind .OR. &
@@ -648,7 +648,7 @@ subroutine read_atms(mype,val_tovs,ithin,isfcalc,&
              idomsfc(1),sfcpct,ts,tsavg,vty,vfr,sty,stp,sm,sn,zz,ff10,sfcr)
      endif
 
-     crit1 = crit1 + rlndsea(isflg) + 10._r_kind*float(iskip) + 0.01_r_kind * abs(zz)
+     crit1 = crit1 + rlndsea(isflg) + 10._r_kind*real(iskip,r_kind) + 0.01_r_kind * abs(zz)
      call checkob(dist1,crit1,itx,iuse)
      if(.not. iuse)cycle ObsLoop
 
@@ -726,7 +726,7 @@ subroutine read_atms(mype,val_tovs,ithin,isfcalc,&
      endif
 
 ! Re-calculate look angle
-     panglr=(start+float(ifov-1)*step)*deg2rad
+     panglr=(start+real(ifov-1,r_kind)*step)*deg2rad
 
 
 !     Load selected observation into data array
diff --git a/src/gsi/read_bufrtovs.f90 b/src/gsi/read_bufrtovs.f90
index a819acd2c3..2fc14b5cdf 100644
--- a/src/gsi/read_bufrtovs.f90
+++ b/src/gsi/read_bufrtovs.f90
@@ -673,7 +673,7 @@ subroutine read_bufrtovs(mype,val_tovs,ithin,isfcalc,&
            terrain = 50._r_kind
            if(llll == 1)terrain = 0.01_r_kind*abs(bfr1bhdr(13))                   
            crit0 = 0.01_r_kind + terrain
-           if (llll >  1 ) crit0 = crit0 + r100 * float(llll)
+           if (llll >  1 ) crit0 = crit0 + r100 * real(llll,r_kind)
            timeinflat=two
            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)
@@ -699,7 +699,7 @@ subroutine read_bufrtovs(mype,val_tovs,ithin,isfcalc,&
            if(hirs .and. ((jsatid == 'n16') .or. (jsatid == 'n17'))) &
               ifovmod=ifovmod+1
 
-           panglr=(start+float(ifovmod-1)*step)*deg2rad
+           panglr=(start+real(ifovmod-1,r_kind)*step)*deg2rad
            lzaest = asin(rato*sin(panglr))
            if( msu .or. hirs2 .or. ssu)then
               lza = lzaest
@@ -821,7 +821,7 @@ subroutine read_bufrtovs(mype,val_tovs,ithin,isfcalc,&
            end do
            if (iskip >= nchanl) cycle read_loop
 !          Map obs to thinning grid
-           crit1 = crit1 + 10._r_kind*float(iskip)
+           crit1 = crit1 + 10._r_kind*real(iskip,r_kind)
            call checkob(dist1,crit1,itx,iuse)
            if(.not. iuse)cycle read_loop
 
diff --git a/src/gsi/read_cris.f90 b/src/gsi/read_cris.f90
index a257480c9a..9843f919d7 100644
--- a/src/gsi/read_cris.f90
+++ b/src/gsi/read_cris.f90
@@ -582,7 +582,7 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,&
 !          Increment nread counter by bufr_nchan    (should be changed to number of channels in satinfo file? (satinfo_nchan))
            nread = nread + satinfo_nchan
            crit0 = 0.01_r_kind
-           if( llll > 1 ) crit0 = crit0 + r100 * float(llll)
+           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)
@@ -600,8 +600,8 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,&
            if( ifor <= 15 )  sat_zenang = -sat_zenang
 
 !          Compute scan angle including sensor twist. 
-           look_angle_est = (start + float((ifor-1))*step) * deg2rad + &
-              fov_dist(ifov) * sin(fov_ang(ifov) - float(ifor-1)*step*deg2rad)
+           look_angle_est = (start + real((ifor-1),r_kind)*step) * deg2rad + &
+              fov_dist(ifov) * sin(fov_ang(ifov) - real(ifor-1,r_kind)*step*deg2rad)
 
            sat_look_angle=asin(rato*sin(sat_zenang*deg2rad))
            if(abs(sat_look_angle)*rad2deg > MAX_SENSOR_ZENITH_ANGLE) then
@@ -763,7 +763,7 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,&
            if(iskip > 0 .and. print_verbose)write(6,*) ' READ_CRIS : iskip > 0 ',iskip
 !          if( iskip >= 10 )cycle read_loop 
 
-           crit1=crit1 + ten*float(iskip)
+           crit1=crit1 + ten*real(iskip,r_kind)
 
 !          Final map obs to grids
            if ( clear ) then 
diff --git a/src/gsi/read_dbz_netcdf.f90 b/src/gsi/read_dbz_netcdf.f90
index 6ea03afaff..ecfb9169c4 100644
--- a/src/gsi/read_dbz_netcdf.f90
+++ b/src/gsi/read_dbz_netcdf.f90
@@ -170,9 +170,9 @@ subroutine read_dbz_mrms_netcdf(nread,ndata,nodata,infile,obstype,lunout,sis,nob
  integer(i_kind) :: num_missing=0,num_nopcp=0, &      !counts 
                     numbadtime=0,num_badtilt=0, &    
                     num_badrange=0,num_m2nopcp=0, &
-		    num_noise=0,num_limmax=0 ,num_limmin=0    
-							     
-        	                                             							     
+                    num_noise=0,num_limmax=0 ,num_limmin=0    
+     
+                                                          
 
 !--General declarations
   integer(i_kind) :: ierror,lunrad,i,j,k,v,na,nb,nelv,nvol, &
@@ -185,7 +185,7 @@ subroutine read_dbz_mrms_netcdf(nread,ndata,nodata,infile,obstype,lunout,sis,nob
   real(r_kind) :: thistiltr,selev0,celev0,thisrange,this_stahgt,thishgt                           
   real(r_kind) :: celev,selev,gamma,thisazimuthr,rlon0, &
                   clat0,slat0,dlat,dlon,thiserr,thislon,thislat, &
-		  rlonloc,rlatloc,rlonglob,rlatglob,timeb,rad_per_meter	  
+                  rlonloc,rlatloc,rlonglob,rlatglob,timeb,rad_per_meter  
   real(r_kind) :: radartwindow
   real(r_kind) :: dbzerr,rmins_an,rmins_ob                                                     
   real(r_kind),allocatable,dimension(:,:):: cdata_all
@@ -203,7 +203,7 @@ subroutine read_dbz_mrms_netcdf(nread,ndata,nodata,infile,obstype,lunout,sis,nob
   !---------SETTINGS FOR FUTURE NAMELIST---------!
   integer(i_kind) :: maxobrange=999000 ! Range (m) *within* which to use observations - obs *outside* this range are not used
   integer(i_kind) :: minobrange=-999 ! Range (m) *outside* of which to use observatons - obs *inside* this range are not used
-  real(r_kind)    :: mintilt=0.0_r_kind   	 ! Only use tilt(elevation) angles (deg) >= this number 
+  real(r_kind)    :: mintilt=0.0_r_kind          ! Only use tilt(elevation) angles (deg) >= this number 
   real(r_kind)    :: maxtilt=20.0_r_kind         ! Do no use tilt(elevation) angles (deg) >= this number
   logical         :: missing_to_nopcp=.false.    ! Set missing observations to 'no precipitation' observations -> dbznoise (See Aksoy et al. 2009, MWR) 
   real(r_kind)    :: dbznoise=2.0_r_kind           ! dBZ obs must be >= dbznoise for assimilation
@@ -241,13 +241,13 @@ subroutine read_dbz_mrms_netcdf(nread,ndata,nodata,infile,obstype,lunout,sis,nob
      if(trim(obstype) == trim(ioctype(i)) .and. abs(icuse(i))== 1) then
         ikx=i 
         radartwindow=ctwind(ikx)*r60         !Time window units converted to minutes 
-	                                     !  (default setting for dbz within convinfo is 0.05 hours)
-	dbzerr=5_r_kind                      !Ob error (dB) to use for radar reflectivity factor
-	exit                                 !Exit loop when finished with initial convinfo fields     
+                                             !  (default setting for dbz within convinfo is 0.05 hours)
+        dbzerr=5_r_kind                      !Ob error (dB) to use for radar reflectivity factor
+        exit                                 !Exit loop when finished with initial convinfo fields     
      else if ( i==nconvtype ) then
         write(6,*) 'READ_dBZ: ERROR - OBSERVATION TYPE IS NOT PRESENT IN CONVINFO OR USE FLAG IS ZERO'
-	write(6,*) 'READ_dBZ: ABORTTING read_dbz.f90 - NO REFLECTIVITY OBS READ!'
-	return 	  
+        write(6,*) 'READ_dBZ: ABORTTING read_dbz.f90 - NO REFLECTIVITY OBS READ!'
+        return   
      endif
   end do     
     
@@ -420,12 +420,12 @@ subroutine read_dbz_mrms_netcdf(nread,ndata,nodata,infile,obstype,lunout,sis,nob
      
         obdate(1)=strct_in_dbz(v,k)%year
         obdate(2)=strct_in_dbz(v,k)%month  
-        obdate(3)=strct_in_dbz(v,k)%day	 
+        obdate(3)=strct_in_dbz(v,k)%day 
         obdate(4)=strct_in_dbz(v,k)%hour   
         obdate(5)=strct_in_dbz(v,k)%minute 
         call w3fs21(obdate,mins_ob)                             !mins_ob -integer number of mins snce 01/01/1978
-	rmins_ob=mins_ob                                        !convert to real number
-	rmins_ob=rmins_ob+(strct_in_dbz(v,k)%second*r60inv)     !convert seconds to minutes and add to ob time
+        rmins_ob=mins_ob                                        !convert to real number
+        rmins_ob=rmins_ob+(strct_in_dbz(v,k)%second*r60inv)     !convert seconds to minutes and add to ob time
  
       !-Comparison is done in units of minutes
       
@@ -439,139 +439,139 @@ subroutine read_dbz_mrms_netcdf(nread,ndata,nodata,infile,obstype,lunout,sis,nob
         if (thistilt <= maxtilt .and. thistilt >= mintilt) then 
      
           gates: do i=1,strct_in_dbz(v,k)%num_gate
-           	  
-   	      thisrange=strct_in_dbz(v,k)%fstgatdis + float(i-1)*strct_in_dbz(v,k)%gateWidth
+             
+             thisrange=strct_in_dbz(v,k)%fstgatdis + real(i-1,r_kind)*strct_in_dbz(v,k)%gateWidth
        
              !-Check to make sure observations are within specified range 
 
-              if (thisrange <= maxobrange .and. thisrange >= minobrange) then	    
-	      azms: do j=1,strct_in_dbz(v,k)%num_beam
-	   
-	           !-Check to see if this is a missing observation 
-		    
-		    nread=nread+1
-		 
+              if (thisrange <= maxobrange .and. thisrange >= minobrange) then    
+                azms: do j=1,strct_in_dbz(v,k)%num_beam
+   
+                    !-Check to see if this is a missing observation 
+    
+                    nread=nread+1
+ 
                     if ( abs(strct_in_dbz(v,k)%field(i,j)) >= 99.0_r_kind ) then
-	               		      
-		      !--Extend no precip observations to missing data fields?
+                     
+                      !--Extend no precip observations to missing data fields?
                       !  May help suppress spurious convection if a problem.
-	  	   
+     
                        if (missing_to_nopcp) then
-	                  strct_in_dbz(v,k)%field(i,j) = dbznoise
-	     	          num_m2nopcp = num_m2nopcp+1
-	               else		  
-			  num_missing=num_missing+1
-	  	          cycle azms                        !No reason to process the ob if it is missing 
+                          strct_in_dbz(v,k)%field(i,j) = dbznoise
+                          num_m2nopcp = num_m2nopcp+1
+                       else  
+                          num_missing=num_missing+1
+                          cycle azms                        !No reason to process the ob if it is missing 
                        end if
-		       		       	         
+                       
                     end if
-		    
-		    
-		    if (l_limmax) then
-		       if ( strct_in_dbz(v,k)%field(i,j) > 60_r_kind ) then
-		          strct_in_dbz(v,k)%field(i,j) = 60_r_kind
-		          num_limmax=num_limmax+1
-		       end if    
-		    end if
-		    if (l_limmin) then
-		       if ( strct_in_dbz(v,k)%field(i,j) < 0_r_kind ) then
-		          strct_in_dbz(v,k)%field(i,j) = 0_r_kind
-		          num_limmin=num_limmin+1
-		       end if    
-		    end if
-		    		    
-		   !--Find observation height using method from read_l2bufr_mod.f90										       
-	         
-		    this_stahgt=strct_in_dbz(v,k)%radhgt
+ 
+    
+                    if (l_limmax) then
+                       if ( strct_in_dbz(v,k)%field(i,j) > 60_r_kind ) then
+                          strct_in_dbz(v,k)%field(i,j) = 60_r_kind
+                          num_limmax=num_limmax+1
+                       end if    
+                    end if
+                    if (l_limmin) then
+                       if ( strct_in_dbz(v,k)%field(i,j) < 0_r_kind ) then
+                          strct_in_dbz(v,k)%field(i,j) = 0_r_kind
+                          num_limmin=num_limmin+1
+                       end if    
+                    end if
+       
+                    !--Find observation height using method from read_l2bufr_mod.f90										       
+         
+                    this_stahgt=strct_in_dbz(v,k)%radhgt
                     aactual=rearth+this_stahgt                    
                     a43=four_thirds*aactual
                     thistiltr=thistilt*deg2rad
                     selev0=sin(thistiltr)
-                    celev0=cos(thistiltr)   		   
-		    b=thisrange*(thisrange+two*aactual*selev0)
+                    celev0=cos(thistiltr)      
+                    b=thisrange*(thisrange+two*aactual*selev0)
                     c=sqrt(aactual*aactual+b)
                     ha=b/(aactual+c)
                     epsh=(thisrange*thisrange-ha*ha)/(r8*aactual)
                     h=ha-epsh
-	            thishgt=this_stahgt+h 
-
-		   !--Find observation location using method from read_l2bufr_mod.f90
-		 
-		   !-Get corrected tilt angle
-	            celev=celev0
-	            selev=selev0
-           	    celev=a43*celev0/(a43+h)
-		    selev=(thisrange*thisrange+h*h+two*a43*h)/(two*thisrange*(a43+h))
-	          
-		    gamma=half*thisrange*(celev0+celev)
-	         
+                    thishgt=this_stahgt+h 
+
+                    !--Find observation location using method from read_l2bufr_mod.f90
+ 
+                    !-Get corrected tilt angle
+                    celev=celev0
+                    selev=selev0
+                    celev=a43*celev0/(a43+h)
+                    selev=(thisrange*thisrange+h*h+two*a43*h)/(two*thisrange*(a43+h))
+          
+                    gamma=half*thisrange*(celev0+celev)
+         
                    !-Get earth lat lon of observation
-	         	 
+          
                     rlon0=deg2rad*strct_in_dbz(v,k)%radlon
-	       	    clat0=cos(deg2rad*strct_in_dbz(v,k)%radlat)
-		    slat0=sin(deg2rad*strct_in_dbz(v,k)%radlat)		  
-		    thisazimuthr=(90.0_r_kind-strct_in_dbz(v,k)%azim(j))*deg2rad   !Storing as 90-azm to
-		                                                                   ! be consistent with 
-										   ! read_l2bufr_mod.f90
-		    rad_per_meter=one/rearth
-		    rlonloc=rad_per_meter*gamma*cos(thisazimuthr)
+                    clat0=cos(deg2rad*strct_in_dbz(v,k)%radlat)
+                    slat0=sin(deg2rad*strct_in_dbz(v,k)%radlat)  
+                    thisazimuthr=(90.0_r_kind-strct_in_dbz(v,k)%azim(j))*deg2rad   !Storing as 90-azm to
+                                                                                   ! be consistent with 
+                                                                                   ! read_l2bufr_mod.f90
+                    rad_per_meter=one/rearth
+                    rlonloc=rad_per_meter*gamma*cos(thisazimuthr)
                     rlatloc=rad_per_meter*gamma*sin(thisazimuthr)
-		                  
-		    call invtllv(rlonloc,rlatloc,rlon0,clat0,slat0,rlonglob,rlatglob)
+                  
+                    call invtllv(rlonloc,rlatloc,rlon0,clat0,slat0,rlonglob,rlatglob)
                  
-		    thislat=rlatglob*rad2deg
+                    thislat=rlatglob*rad2deg
                     thislon=rlonglob*rad2deg 
   
                    !-Check format of longitude and correct if necessary
                  
-		    if(thislon>=r360) thislon=thislon-r360
+                    if(thislon>=r360) thislon=thislon-r360
                     if(thislon<zero ) thislon=thislon+r360
                  
-		   !-Convert back to radians                 
-         		       		       
-		    thislat = thislat*deg2rad
+                   !-Convert back to radians                 
+                      
+                    thislat = thislat*deg2rad
                     thislon = thislon*deg2rad
-                 		 
-		    !find grid relative lat lon locations of earth lat lon
+                  
+                   !find grid relative lat lon locations of earth lat lon
                  
-		    call tll2xy(thislon,thislat,dlon,dlat,outside)
+                    call tll2xy(thislon,thislat,dlon,dlat,outside)
                     if (outside) cycle azms             !If observation is outside the domain
-		                                        ! then cycle, but don't increase range right away.
-							! Domain could be rectangular, so ob may be out of
-						        ! range at one end, but not the other.		     					                   		   		   
-		   thiserr=dbzerr
+                                                        ! then cycle, but don't increase range right away.
+                                                        ! Domain could be rectangular, so ob may be out of
+                                                        ! range at one end, but not the other.		     					                   		   		   
+                    thiserr=dbzerr
                 
-		    	     
-		    ndata  = min(ndata+1,maxobs)     
-		    nodata = min(nodata+1,maxobs)  !number of obs not used (no meaning here)
-		    
-		    this_staid=strct_in_dbz(v,k)%radid      !Via equivalence in declaration, value is propagated
-		    					    !  to rstation_id used below.
-		    
-		    cdata_all(1,ndata) = thiserr 		      ! reflectivity obs error (dB) - inflated/adjusted
-		    cdata_all(2,ndata) = dlon			      ! grid relative longitude
-		    cdata_all(3,ndata) = dlat			      ! grid relative latitude
-		    cdata_all(4,ndata) = thishgt		      ! obs absolute height (m)
-		    cdata_all(5,ndata) = strct_in_dbz(v,k)%field(i,j) ! radar reflectivity factor 
-		    cdata_all(6,ndata) = thisazimuthr		      ! 90deg-azimuth angle (radians)
-		    cdata_all(7,ndata) = timeb*r60inv		      ! obs time (analyis relative hour)
-		    cdata_all(8,ndata) = ikx			      ! type	 double check with the convinfo txt	   
-		    cdata_all(9,ndata) = thistiltr		      ! tilt angle (radians)
-		    cdata_all(10,ndata)= this_stahgt		      ! station elevation (m)
-		    cdata_all(11,ndata)= rstation_id		      ! station id
-		    cdata_all(12,ndata)= icuse(ikx)		      ! usage parameter
-		    cdata_all(13,ndata)= thislon*rad2deg	      ! earth relative longitude (degrees)
-		    cdata_all(14,ndata)= thislat*rad2deg	      ! earth relative latitude (degrees)
-		    cdata_all(15,ndata)= thisrange		      ! range from radar in m 
-		    cdata_all(16,ndata)= dbzerr                       ! orginal error from convinfo file
-		    cdata_all(17,ndata)= dbznoise                     ! noise threshold for reflectivity (dBZ)
+         
+                    ndata  = min(ndata+1,maxobs)     
+                    nodata = min(nodata+1,maxobs)  !number of obs not used (no meaning here)
+    
+                    this_staid=strct_in_dbz(v,k)%radid      !Via equivalence in declaration, value is propagated
+                                                            !  to rstation_id used below.
+   
+                    cdata_all(1,ndata) = thiserr                      ! reflectivity obs error (dB) - inflated/adjusted
+                    cdata_all(2,ndata) = dlon                         ! grid relative longitude
+                    cdata_all(3,ndata) = dlat                         ! grid relative latitude
+                    cdata_all(4,ndata) = thishgt                      ! obs absolute height (m)
+                    cdata_all(5,ndata) = strct_in_dbz(v,k)%field(i,j) ! radar reflectivity factor 
+                    cdata_all(6,ndata) = thisazimuthr                 ! 90deg-azimuth angle (radians)
+                    cdata_all(7,ndata) = timeb*r60inv                 ! obs time (analyis relative hour)
+                    cdata_all(8,ndata) = ikx                          ! type double check with the convinfo txt	   
+                    cdata_all(9,ndata) = thistiltr                    ! tilt angle (radians)
+                    cdata_all(10,ndata)= this_stahgt                  ! station elevation (m)
+                    cdata_all(11,ndata)= rstation_id                  ! station id
+                    cdata_all(12,ndata)= icuse(ikx)                   ! usage parameter
+                    cdata_all(13,ndata)= thislon*rad2deg              ! earth relative longitude (degrees)
+                    cdata_all(14,ndata)= thislat*rad2deg              ! earth relative latitude (degrees)
+                    cdata_all(15,ndata)= thisrange                    ! range from radar in m 
+                    cdata_all(16,ndata)= dbzerr                       ! orginal error from convinfo file
+                    cdata_all(17,ndata)= dbznoise                     ! noise threshold for reflectivity (dBZ)
 
                   end do azms  !j
               else
-	         num_badrange=num_badrange+1      !If outside acceptable range, increment
-	      end if   !Range check	
-		
-	   end do gates    !i
+                  num_badrange=num_badrange+1      !If outside acceptable range, increment
+              end if   !Range check	
+
+           end do gates    !i
      
         else
            num_badtilt=num_badtilt+1           !If outside acceptable tilts, increment
@@ -794,9 +794,9 @@ subroutine read_dbz_mrms_sparse_netcdf(nread,ndata,nodata,infile,obstype,lunout,
  integer(i_kind) :: num_missing=0,num_nopcp=0, &      !counts 
                     numbadtime=0,num_badtilt=0, &    
                     num_badrange=0,num_m2nopcp=0, &
-		    num_noise=0,num_limmax=0 ,num_limmin=0   
-							     
-        	                                             							     
+                    num_noise=0,num_limmax=0 ,num_limmin=0   
+    
+                                                          
 
 !--General declarations
   integer(i_kind) :: ierror,lunrad,i,j,k,v,na,nb,nelv,nvol, &
@@ -809,7 +809,7 @@ subroutine read_dbz_mrms_sparse_netcdf(nread,ndata,nodata,infile,obstype,lunout,
   real(r_kind) :: thistiltr,selev0,celev0,thisrange,this_stahgt,thishgt                           
   real(r_kind) :: celev,selev,gamma,thisazimuthr,rlon0, &
                   clat0,slat0,dlat,dlon,thiserr,thislon,thislat, &
-		  rlonloc,rlatloc,rlonglob,rlatglob,timeb,rad_per_meter	  
+                  rlonloc,rlatloc,rlonglob,rlatglob,timeb,rad_per_meter  
   real(r_kind) :: radartwindow
   real(r_kind) :: dbzerr,rmins_an,rmins_ob                                                     
   real(r_kind),allocatable,dimension(:,:):: cdata_all
@@ -827,7 +827,7 @@ subroutine read_dbz_mrms_sparse_netcdf(nread,ndata,nodata,infile,obstype,lunout,
   !---------SETTINGS FOR FUTURE NAMELIST---------!
   integer(i_kind) :: maxobrange=99900000 ! Range (m) *within* which to use observations - obs *outside* this range are not used
   integer(i_kind) :: minobrange=-999 ! Range (m) *outside* of which to use observatons - obs *inside* this range are not used
-  real(r_kind)    :: mintilt=0.0_r_kind   	 ! Only use tilt(elevation) angles (deg) >= this number 
+  real(r_kind)    :: mintilt=0.0_r_kind          ! Only use tilt(elevation) angles (deg) >= this number 
   real(r_kind)    :: maxtilt=20.0_r_kind         ! Do no use tilt(elevation) angles (deg) >= this number
   logical         :: missing_to_nopcp=.false.    ! Set missing observations to 'no precipitation' observations -> dbznoise (See Aksoy et al. 2009, MWR) 
   real(r_kind)    :: dbznoise=2_r_kind           ! dBZ obs must be >= dbznoise for assimilation
@@ -863,13 +863,13 @@ subroutine read_dbz_mrms_sparse_netcdf(nread,ndata,nodata,infile,obstype,lunout,
      if(trim(obstype) == trim(ioctype(i)) .and. abs(icuse(i))== 1) then
         ikx=i 
         radartwindow=ctwind(ikx)*r60         !Time window units converted to minutes 
-	                                     !  (default setting for dbz within convinfo is 0.05 hours)
-	dbzerr=5_r_kind                      !Ob error (dB) to use for radar reflectivity factor
-	exit                                 !Exit loop when finished with initial convinfo fields     
+                                             !  (default setting for dbz within convinfo is 0.05 hours)
+        dbzerr=5_r_kind                      !Ob error (dB) to use for radar reflectivity factor
+        exit                                 !Exit loop when finished with initial convinfo fields     
      else if ( i==nconvtype ) then
         write(6,*) 'READ_dBZ: ERROR - OBSERVATION TYPE IS NOT PRESENT IN CONVINFO OR USE FLAG IS ZERO'
-	write(6,*) 'READ_dBZ: ABORTTING read_dbz.f90 - NO REFLECTIVITY OBS READ!'
-	return 	  
+        write(6,*) 'READ_dBZ: ABORTTING read_dbz.f90 - NO REFLECTIVITY OBS READ!'
+        return   
      endif
   end do     
     
@@ -1017,20 +1017,20 @@ subroutine read_dbz_mrms_sparse_netcdf(nread,ndata,nodata,infile,obstype,lunout,
 
 ! transform the read-in ob to the intermidate  obs variables( radar obs  to be used in GSI
 
-         strct_in_dbz(v,k)%radid=radarsite_nc
-         strct_in_dbz(v,k)%vcpnum=vcp_nc
-         strct_in_dbz(v,k)%year=iyear  !  to be defind from infile name              
+        strct_in_dbz(v,k)%radid=radarsite_nc
+        strct_in_dbz(v,k)%vcpnum=vcp_nc
+        strct_in_dbz(v,k)%year=iyear  !  to be defind from infile name              
         strct_in_dbz(v,k)%month=imon                
         strct_in_dbz(v,k)%day=iday                  
         strct_in_dbz(v,k)%hour=ihour                 
         strct_in_dbz(v,k)%minute=imin               
         strct_in_dbz(v,k)%second=isec
-         strct_in_dbz(v,k)%radlat=lat_nc
+        strct_in_dbz(v,k)%radlat=lat_nc
         strct_in_dbz(v,k)%radlon=lon_nc   
         strct_in_dbz(v,k)%radhgt=height_nc
-         strct_in_dbz(v,k)%fstgatdis =firstgate_nc      
+        strct_in_dbz(v,k)%fstgatdis =firstgate_nc      
         strct_in_dbz(v,k)%gateWidth=gatewidth_nc(1) ! always the same ??)
-         strct_in_dbz(v,k)%elev_angle=elev_nc
+        strct_in_dbz(v,k)%elev_angle=elev_nc
         strct_in_dbz(v,k)%num_beam=numazim_nc          
         strct_in_dbz(v,k)%num_gate=numgate_nc
         na=strct_in_dbz(v,k)%num_beam
@@ -1065,8 +1065,8 @@ subroutine read_dbz_mrms_sparse_netcdf(nread,ndata,nodata,infile,obstype,lunout,
         obdate(5)=strct_in_dbz(v,k)%minute 
         call w3fs21(obdate,mins_ob)                             !mins_ob -integer number of mins snce 01/01/1978
        
-	rmins_ob=mins_ob                                        !convert to real number
-	rmins_ob=rmins_ob+(strct_in_dbz(v,k)%second*r60inv)     !convert seconds to minutes and add to ob time
+        rmins_ob=mins_ob                                        !convert to real number
+        rmins_ob=rmins_ob+(strct_in_dbz(v,k)%second*r60inv)     !convert seconds to minutes and add to ob time
  
       !-Comparison is done in units of minutes
       
@@ -1088,139 +1088,139 @@ subroutine read_dbz_mrms_sparse_netcdf(nread,ndata,nodata,infile,obstype,lunout,
           pixel: do ipix=1,real_numpixel
               j=pixel_x_nc(ipix)+1
               i=pixel_y_nc(ipix)+1
-           	  
-   	      thisrange=strct_in_dbz(v,k)%fstgatdis + float(i-1)*strct_in_dbz(v,k)%gateWidth
+             
+              thisrange=strct_in_dbz(v,k)%fstgatdis + real(i-1,r_kind)*strct_in_dbz(v,k)%gateWidth
        
              !-Check to make sure observations are within specified range 
 
-              if (thisrange <= maxobrange .and. thisrange >= minobrange) then	    
-	  	   
-		    
-		    nread=nread+1
+              if (thisrange <= maxobrange .and. thisrange >= minobrange) then    
+     
+    
+                    nread=nread+1
                     if ( abs(obdata_pixel_nc(ipix)) >= 999.0_r_kind ) then
-	               		      
-		      !--Extend no precip observations to missing data fields?
+
+                      !--Extend no precip observations to missing data fields?
                       !  May help suppress spurious convection if a problem.
-	  	   
+     
                        if (missing_to_nopcp) then
-	                  obdata_pixel_nc(ipix) = dbznoise
-	     	          num_m2nopcp = num_m2nopcp+1
-	               else		  
-			  num_missing=num_missing+1
-	  	          cycle pixel                       !No reason to process the ob if it is missing 
+                          obdata_pixel_nc(ipix) = dbznoise
+                          num_m2nopcp = num_m2nopcp+1
+                       else  
+                          num_missing=num_missing+1
+                          cycle pixel                       !No reason to process the ob if it is missing 
                        end if
-		       		       	         
+
+                    end if
+
+
+                    if (l_limmax) then
+                       if ( obdata_pixel_nc(ipix) > 60_r_kind ) then
+                          obdata_pixel_nc(ipix) = 60_r_kind
+                          num_limmax=num_limmax+1
+                       end if    
                     end if
-		    
-		    
-		    if (l_limmax) then
-		       if ( obdata_pixel_nc(ipix) > 60_r_kind ) then
-		          obdata_pixel_nc(ipix) = 60_r_kind
-		          num_limmax=num_limmax+1
-		       end if    
-		    end if
-		    if (l_limmin) then
-		       if ( obdata_pixel_nc(ipix) < 0_r_kind ) then
-		          obdata_pixel_nc(ipix) = 0_r_kind
-		          num_limmin=num_limmin+1
-		       end if    
-		    end if
-		    		    
-		   !-Special treatment for no-precip obs?	       			       
-
-                                     			
-		   !--Find observation height using method from read_l2bufr_mod.f90										       
-	         
-		    this_stahgt=strct_in_dbz(v,k)%radhgt
+                    if (l_limmin) then
+                       if ( obdata_pixel_nc(ipix) < 0_r_kind ) then
+                          obdata_pixel_nc(ipix) = 0_r_kind
+                          num_limmin=num_limmin+1
+                       end if    
+                    end if
+
+                   !-Special treatment for no-precip obs?	       			       
+
+
+                   !--Find observation height using method from read_l2bufr_mod.f90										       
+
+                    this_stahgt=strct_in_dbz(v,k)%radhgt
                     aactual=rearth+this_stahgt                    
                     a43=four_thirds*aactual
                     thistiltr=thistilt*deg2rad
                     selev0=sin(thistiltr)
-                    celev0=cos(thistiltr)   		   
-		    b=thisrange*(thisrange+two*aactual*selev0)
+                    celev0=cos(thistiltr)      
+                    b=thisrange*(thisrange+two*aactual*selev0)
                     c=sqrt(aactual*aactual+b)
                     ha=b/(aactual+c)
                     epsh=(thisrange*thisrange-ha*ha)/(r8*aactual)
                     h=ha-epsh
-	            thishgt=this_stahgt+h 
-
-		   !--Find observation location using method from read_l2bufr_mod.f90
-		 
-		   !-Get corrected tilt angle
-	            celev=celev0
-	            selev=selev0
-           	    celev=a43*celev0/(a43+h)
-		    selev=(thisrange*thisrange+h*h+two*a43*h)/(two*thisrange*(a43+h))
-	          
-		    gamma=half*thisrange*(celev0+celev)
-	         
+                    thishgt=this_stahgt+h 
+
+                   !--Find observation location using method from read_l2bufr_mod.f90
+ 
+                   !-Get corrected tilt angle
+                    celev=celev0
+                    selev=selev0
+                    celev=a43*celev0/(a43+h)
+                    selev=(thisrange*thisrange+h*h+two*a43*h)/(two*thisrange*(a43+h))
+          
+                    gamma=half*thisrange*(celev0+celev)
+         
                    !-Get earth lat lon of observation
-	         	 
+          
                     rlon0=deg2rad*strct_in_dbz(v,k)%radlon
-	       	    clat0=cos(deg2rad*strct_in_dbz(v,k)%radlat)
-		    slat0=sin(deg2rad*strct_in_dbz(v,k)%radlat)		  
-		    thisazimuthr=(90.0_r_kind-strct_in_dbz(v,k)%azim(j))*deg2rad   !Storing as 90-azm to
-		                                                                   ! be consistent with 
-										   ! read_l2bufr_mod.f90
-		    rad_per_meter=one/rearth
-		    rlonloc=rad_per_meter*gamma*cos(thisazimuthr)
+                    clat0=cos(deg2rad*strct_in_dbz(v,k)%radlat)
+                    slat0=sin(deg2rad*strct_in_dbz(v,k)%radlat)  
+                    thisazimuthr=(90.0_r_kind-strct_in_dbz(v,k)%azim(j))*deg2rad   !Storing as 90-azm to
+                                                                                   ! be consistent with 
+                                                                                   ! read_l2bufr_mod.f90
+                    rad_per_meter=one/rearth
+                    rlonloc=rad_per_meter*gamma*cos(thisazimuthr)
                     rlatloc=rad_per_meter*gamma*sin(thisazimuthr)
-		                  
-		    call invtllv(rlonloc,rlatloc,rlon0,clat0,slat0,rlonglob,rlatglob)
+                  
+                    call invtllv(rlonloc,rlatloc,rlon0,clat0,slat0,rlonglob,rlatglob)
                  
-		    thislat=rlatglob*rad2deg
+                    thislat=rlatglob*rad2deg
                     thislon=rlonglob*rad2deg 
   
                    !-Check format of longitude and correct if necessary
                  
-		    if(thislon>=r360) thislon=thislon-r360
+                    if(thislon>=r360) thislon=thislon-r360
                     if(thislon<zero ) thislon=thislon+r360
                  
-		   !-Convert back to radians                 
-         		       		       
-		    thislat = thislat*deg2rad
+                   !-Convert back to radians                 
+                      
+                    thislat = thislat*deg2rad
                     thislon = thislon*deg2rad
-                 		 
-		    !find grid relative lat lon locations of earth lat lon
                  
-		    call tll2xy(thislon,thislat,dlon,dlat,outside)
+                   !find grid relative lat lon locations of earth lat lon
+                 
+                    call tll2xy(thislon,thislat,dlon,dlat,outside)
                     if (outside) cycle pixel             !If observation is outside the domain
-		   thiserr=dbzerr
+                    thiserr=dbzerr
                 
-		      
-		      
-		      
+
+
+
 
                    !-Load good data into output array
-		    	     
-		    ndata  = min(ndata+1,maxobs)     
-		    nodata = min(nodata+1,maxobs)  !number of obs not used (no meaning here)
-		    
-		    this_staid=strct_in_dbz(v,k)%radid      !Via equivalence in declaration, value is propagated
-		    					    !  to rstation_id used below.
-		    
-		    cdata_all(1,ndata) = thiserr 		      ! reflectivity obs error (dB) - inflated/adjusted
-		    cdata_all(2,ndata) = dlon			      ! grid relative longitude
-		    cdata_all(3,ndata) = dlat			      ! grid relative latitude
-		    cdata_all(4,ndata) = thishgt		      ! obs absolute height (m)
-		    cdata_all(5,ndata) = obdata_pixel_nc(ipix)  !strct_in_dbz(v,k)%field(i,j) ! radar reflectivity factor 
-		    cdata_all(6,ndata) = thisazimuthr		      ! 90deg-azimuth angle (radians)
-		    cdata_all(7,ndata) = timeb*r60inv		      ! obs time (analyis relative hour)
-		    cdata_all(8,ndata) = ikx			      ! type		   
-		    cdata_all(9,ndata) = thistiltr		      ! tilt angle (radians)
-		    cdata_all(10,ndata)= this_stahgt		      ! station elevation (m)
-		    cdata_all(11,ndata)= rstation_id		      ! station id
-		    cdata_all(12,ndata)= icuse(ikx)		      ! usage parameter
-		    cdata_all(13,ndata)= thislon*rad2deg	      ! earth relative longitude (degrees)
-		    cdata_all(14,ndata)= thislat*rad2deg	      ! earth relative latitude (degrees)
-		    cdata_all(15,ndata)= thisrange		      ! range from radar in m 
-		    cdata_all(16,ndata)= dbzerr                       ! orginal error from convinfo file
-		    cdata_all(17,ndata)= dbznoise                     ! noise threshold for reflectivity (dBZ)
+         
+                    ndata  = min(ndata+1,maxobs)     
+                    nodata = min(nodata+1,maxobs)  !number of obs not used (no meaning here)
+    
+                    this_staid=strct_in_dbz(v,k)%radid      !Via equivalence in declaration, value is propagated
+                                                            !  to rstation_id used below.
+
+                    cdata_all(1,ndata) = thiserr                      ! reflectivity obs error (dB) - inflated/adjusted
+                    cdata_all(2,ndata) = dlon                         ! grid relative longitude
+                    cdata_all(3,ndata) = dlat                         ! grid relative latitude
+                    cdata_all(4,ndata) = thishgt                      ! obs absolute height (m)
+                    cdata_all(5,ndata) = obdata_pixel_nc(ipix)  !strct_in_dbz(v,k)%field(i,j) ! radar reflectivity factor 
+                    cdata_all(6,ndata) = thisazimuthr                 ! 90deg-azimuth angle (radians)
+                    cdata_all(7,ndata) = timeb*r60inv                 ! obs time (analyis relative hour)
+                    cdata_all(8,ndata) = ikx                          ! type		   
+                    cdata_all(9,ndata) = thistiltr                    ! tilt angle (radians)
+                    cdata_all(10,ndata)= this_stahgt                  ! station elevation (m)
+                    cdata_all(11,ndata)= rstation_id                  ! station id
+                    cdata_all(12,ndata)= icuse(ikx)                   ! usage parameter
+                    cdata_all(13,ndata)= thislon*rad2deg              ! earth relative longitude (degrees)
+                    cdata_all(14,ndata)= thislat*rad2deg              ! earth relative latitude (degrees)
+                    cdata_all(15,ndata)= thisrange                    ! range from radar in m 
+                    cdata_all(16,ndata)= dbzerr                       ! orginal error from convinfo file
+                    cdata_all(17,ndata)= dbznoise                     ! noise threshold for reflectivity (dBZ)
 
               else
-	         num_badrange=num_badrange+1      !If outside acceptable range, increment
-	      end if   !Range check	
-		
+                   num_badrange=num_badrange+1      !If outside acceptable range, increment
+              end if   !Range check	
+
      
               end do pixel    !i
         else
diff --git a/src/gsi/read_files.f90 b/src/gsi/read_files.f90
index dadcbff3e5..60282bfbc6 100644
--- a/src/gsi/read_files.f90
+++ b/src/gsi/read_files.f90
@@ -69,7 +69,8 @@ subroutine read_files(mype)
 !                                      nfsecondn  FCST Secs (i_kind) numerator
 !                                      nfsecondd  FCST Secs (i_kind) denominator
 !
-!       %fhour = float(nfhour) + float(nfminute)/r60 + float(nfsecondn)/float(nfsecondd)/r3600
+!       %fhour = real(nfhour,r_kind) + real(nfminute,r_kind)/r60 + &
+!                real(nfsecondn,r_kind)/real(nfsecondd,r_kind)/r3600
 
 ! attributes:
 !   language: f90
@@ -264,7 +265,7 @@ subroutine read_files(mype)
                 write(6,*)'READ_FILES: ***WARNING*** problem reading atm file ',trim(filename),iret 
            idate6 = get_idate_from_time_units(atmges) 
            call read_vardata(atmges, 'time', fhour)
-           hourg4 = float(nint(fhour(1))) ! going to make this nearest integer for now
+           hourg4 = real(nint(fhour(1)),r_kind) ! going to make this nearest integer for now
            idateg(1) = idate6(4)
            idateg(2) = idate6(2)
            idateg(3) = idate6(3)
@@ -288,7 +289,8 @@ subroutine read_files(mype)
               call stop2(80)
            endif
 
-           hourg4 = float(nfhour) + float(nfminute)/r60 + float(nfsecondn)/float(nfsecondd)/r3600
+           hourg4 = real(nfhour,r_kind) + real(nfminute,r_kind)/r60 + &
+                    real(nfsecondn,r_kind)/real(nfsecondd,r_kind)/r3600
            idateg(1) = idate(4)  !hour
            idateg(2) = idate(2)  !month
            idateg(3) = idate(3)  !day
@@ -342,7 +344,7 @@ subroutine read_files(mype)
            ncdim = get_dim(sfcges, 'grid_yt'); sfc_head%latb = ncdim%len
            idate6 = get_idate_from_time_units(sfcges) 
            call read_vardata(sfcges, 'time', fhour)
-           hourg4 = float(nint(fhour(1))) ! going to make this nearest integer for now
+           hourg4 = real(nint(fhour(1)),r_kind) ! going to make this nearest integer for now
            idateg(1) = idate6(4)
            idateg(2) = idate6(2)
            idateg(3) = idate6(3)
@@ -381,7 +383,8 @@ subroutine read_files(mype)
                  nfhour, nfminute, nfsecondn, nfsecondd
               call stop2(80)
            endif
-           hourg4   = float(nfhour) + float(nfminute)/r60 + float(nfsecondn)/float(nfsecondd)/r3600
+           hourg4   = real(nfhour,r_kind) + real(nfminute,r_kind)/r60 + &
+                      real(nfsecondn,r_kind)/real(nfsecondd,r_kind)/r3600
            idateg(1) = idate(4)  !hour
            idateg(2) = idate(2)  !month
            idateg(3) = idate(3)  !day
@@ -478,7 +481,8 @@ subroutine read_files(mype)
                          nfhour, nfminute, nfsecondn, nfsecondd
                     call stop2(80)
                  endif
-                 hourg4   = float(nfhour) + float(nfminute)/r60 + float(nfsecondn)/float(nfsecondd)/r3600
+                 hourg4   = real(nfhour,r_kind) + real(nfminute,r_kind)/r60 + &
+                            real(nfsecondn,r_kind)/real(nfsecondd,r_kind)/r3600
                  idateg(1) = idate(4)  !hour
                  idateg(2) = idate(2)  !month
                  idateg(3) = idate(3)  !day
@@ -540,7 +544,8 @@ subroutine read_files(mype)
                       nfhour, nfminute, nfsecondn, nfsecondd
                  call stop2(80)
               endif
-              hourg4   = float(nfhour) + float(nfminute)/r60 + float(nfsecondn)/float(nfsecondd)/r3600
+              hourg4   = real(nfhour,r_kind) + real(nfminute,r_kind)/r60 + &
+                         real(nfsecondn,r_kind)/real(nfsecondd,r_kind)/r3600
               idateg(1) = idate(4)  !hour
               idateg(2) = idate(2)  !month
               idateg(3) = idate(3)  !day
diff --git a/src/gsi/read_gfs_ozone_for_regional.f90 b/src/gsi/read_gfs_ozone_for_regional.f90
index 5b6acce30b..e018d19cba 100644
--- a/src/gsi/read_gfs_ozone_for_regional.f90
+++ b/src/gsi/read_gfs_ozone_for_regional.f90
@@ -291,8 +291,8 @@ subroutine read_gfs_ozone_for_regional
      if (nframe /= 0) call error_msg(trim(my_name),trim(filename),'nframe', &
                                      'getfilehead',istop,nframe)
 
-     fhour = float(nfhour) + float(nfminute)/r60 + &
-             float(nfsecondn)/float(nfsecondd)/r3600
+     fhour = real(nfhour,r_kind) + real(nfminute,r_kind)/r60 + &
+             real(nfsecondn,r_kind)/real(nfsecondd,r_kind)/r3600
 
      nsig_gfs=levs
      jcap_org=njcap
diff --git a/src/gsi/read_gmi.f90 b/src/gsi/read_gmi.f90
index f59529662a..f54263b129 100644
--- a/src/gsi/read_gmi.f90
+++ b/src/gsi/read_gmi.f90
@@ -678,7 +678,7 @@ subroutine read_gmi(mype,val_gmi,ithin,rmesh,jsatid,gstime,&
        if(.not. regional .and. dist1 > 0.75_r_kind) cycle obsloop
     endif
 
-    crit1 = crit1 + 10._r_kind * float(iskip)
+    crit1 = crit1 + 10._r_kind * real(iskip,r_kind)
     call checkob(dist1,crit1,itx,iuse)
     if(.not. iuse) then
        cycle obsloop
@@ -819,7 +819,7 @@ subroutine read_gmi(mype,val_gmi,ithin,rmesh,jsatid,gstime,&
          if(pos_max==0) then
            j2=1
          else
-           j2=nint(float(pos_statis(i))/pos_max)
+           j2=nint(real(pos_statis(i),r_kind)/pos_max)
            j2=max(1,j2)
          endif
          do j=1,pos_statis(i),j2
@@ -835,7 +835,7 @@ subroutine read_gmi(mype,val_gmi,ithin,rmesh,jsatid,gstime,&
          if(pos_max==0) then
            j2=1
          else
-           j2=nint(float(pos_statis(i))/pos_max)
+           j2=nint(real(pos_statis(i),r_kind)/pos_max)
            j2=max(1,j2)
          endif
          do j=1,pos_statis(i),j2
diff --git a/src/gsi/read_goesglm.f90 b/src/gsi/read_goesglm.f90
index bf8639c72d..f087430092 100644
--- a/src/gsi/read_goesglm.f90
+++ b/src/gsi/read_goesglm.f90
@@ -224,7 +224,7 @@ subroutine read_goesglm(nread,ndata,nodata,infile,obstype,lunout,twindin,sis)
             call w3fs21(idate5,minan)    !  analysis ref time in seconds relative to historic date
            
 !           Add obs reference time, then subtract analysis time to get obs time relative to analysis
-            time_correction=float(minobs-minan)*r60inv
+            time_correction=real(minobs-minan,r_kind)*r60inv
         else
             time_correction=zero
         end if
@@ -512,13 +512,13 @@ subroutine convert_to_flash_rate   &
 
      end do  !! do iobs=2,ndata_strike
 
-     darea=darea_sum/float(ndata_strike)
+     darea=darea_sum/real(ndata_strike,r_kind)
   else   !! ndata_strike=0
      darea=zero
 
   end if  !! if(ndata_strike>0) then
 
-  dtime=float(nhr_assimilation)
+  dtime=real(nhr_assimilation,r_kind)
 
   ! Regional
 
@@ -574,8 +574,8 @@ subroutine convert_to_flash_rate   &
 
 !!   find lightning strikes near the (ii0,jj0) point
 
-     xbound=float(ii0)
-     ybound=float(jj0)
+     xbound=real(ii0,r_kind)
+     ybound=real(jj0,r_kind)
 
      xflag=(xx>xbound) .AND. (xx<xbound+1.)
      yflag=(yy>ybound) .AND. (yy<ybound+1.)
@@ -606,10 +606,10 @@ subroutine convert_to_flash_rate   &
 
   do index=1,ngridh
      if (lcount(index)>0) then
-        glon_central(index)=glon_central(index)/float(lcount(index))
-        glat_central(index)=glat_central(index)/float(lcount(index))
-        lon_central(index)= lon_central(index)/float(lcount(index))
-        lat_central(index)= lat_central(index)/float(lcount(index))
+        glon_central(index)=glon_central(index)/real(lcount(index),r_kind)
+        glat_central(index)=glat_central(index)/real(lcount(index),r_kind)
+        lon_central(index)= lon_central(index)/real(lcount(index),r_kind)
+        lat_central(index)= lat_central(index)/real(lcount(index),r_kind)
      endif  !! if(lcount(index)>0) then
   enddo  !! do index=1,ngridh
 
@@ -653,7 +653,7 @@ subroutine convert_to_flash_rate   &
         cdata_flash_h( 3,icount)=glat_central(index)
 
         if (darea>0._r_kind) then
-           cdata_flash_h( 4,icount)=float(lcount(index))/(darea*dtime)
+           cdata_flash_h( 4,icount)=real(lcount(index),r_kind)/(darea*dtime)
         else 
            cdata_flash_h( 4,icount)=0. 
         end if
@@ -727,22 +727,22 @@ subroutine convert_time (date_old,date_new,nmax)
      jdd=INT(0.0001_r_kind*xdate(i))
      idd=INT(xdate(i))-jdd*10000
 
-     ysumidd=float(idd)
-     dd=float(INT(0.01_r_kind*ysumidd))
+     ysumidd=real(idd,r_kind)
+     dd=real(INT(0.01_r_kind*ysumidd),r_kind)
      hh=ysumidd-dd*100._r_kind
 
      sumidd=sumidd+dd*24._r_kind+hh
 
   enddo  !! do i=1,nmax
 
-  xsumidd=float(sumidd)/nmax
-  ysumidd=float(INT(xsumidd))
+  xsumidd=real(sumidd,r_kind)/nmax
+  ysumidd=real(INT(xsumidd),r_kind)
 
   kdd=INT(xsumidd/24._r_kind)
-  xdd=float(kdd)
-  xhh=ysumidd-float(kdd)*24._r_kind
+  xdd=real(kdd,r_kind)
+  xhh=ysumidd-real(kdd,r_kind)*24._r_kind
 
-  ydate=float(jdd)*10000._r_kind+xdd*100._r_kind+xhh+xccyy
+  ydate=real(jdd,r_kind)*10000._r_kind+xdd*100._r_kind+xhh+xccyy
       
   date_old=ydate 
       
diff --git a/src/gsi/read_goesndr.f90 b/src/gsi/read_goesndr.f90
index 86fc1f0a5c..7c55b6ab4c 100644
--- a/src/gsi/read_goesndr.f90
+++ b/src/gsi/read_goesndr.f90
@@ -382,7 +382,7 @@ subroutine read_goesndr(mype,val_goes,ithin,rmesh,jsatid,infile,&
 
         nread=nread+nchanl
         crit0=0.01_r_kind
-        if(ifov < mfov .and. ifov > 0) crit0 = crit0+two*float(mfov-ifov)
+        if(ifov < mfov .and. ifov > 0) crit0 = crit0+two*real(mfov-ifov,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)
diff --git a/src/gsi/read_iasi.f90 b/src/gsi/read_iasi.f90
index 4e688cd36e..9ab1d5446f 100644
--- a/src/gsi/read_iasi.f90
+++ b/src/gsi/read_iasi.f90
@@ -577,7 +577,7 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,&
            nread = nread + satinfo_nchan
 
            crit0 = 0.01_r_kind
-           if( llll > 1 ) crit0 = crit0 + r100 * float(llll)
+           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)
@@ -598,7 +598,7 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,&
 !          Compare IASI satellite scan angle and zenith angle
            piece = -step_adjust
            if ( mod(ifovn,2) == 1) piece = step_adjust
-           lza = ((start + float((ifov-1)/4)*step) + piece)*deg2rad
+           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
@@ -748,7 +748,7 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,&
               cycle read_loop 
            end if
 
-!          crit1=crit1 + ten*float(iskip)
+!          crit1=crit1 + ten*real(iskip,r_kind)
 
 !          If the surface channel exists (~960.0 cm-1) and the AVHRR cloud information is missing, use an
 !          estimate of the surface temperature to determine if the profile may be clear.
diff --git a/src/gsi/read_lidar.f90 b/src/gsi/read_lidar.f90
index 6d74de0802..ad5b27b784 100644
--- a/src/gsi/read_lidar.f90
+++ b/src/gsi/read_lidar.f90
@@ -172,7 +172,7 @@ subroutine read_lidar(nread,ndata,nodata,infile,obstype,lunout,twind,sis,nobs)
 
 !     add obs reference time, then subtract analysis time to get obs time relative to analysis
 
-     time_correction=float(minobs-minan)*r60inv
+     time_correction=real(minobs-minan,r_kind)*r60inv
 
   else
      time_correction=zero
diff --git a/src/gsi/read_nsstbufr.f90 b/src/gsi/read_nsstbufr.f90
index adfdee4f13..d7a3472dd0 100644
--- a/src/gsi/read_nsstbufr.f90
+++ b/src/gsi/read_nsstbufr.f90
@@ -367,7 +367,7 @@ subroutine read_nsstbufr(nread,ndata,nodata,gstime,infile,obstype,lunout, &
            if ( rsc > 60.0_r_kind .or. rsc < zero ) rsc = zero   !second in real
 
            call w3fs21(idate5,nmind)
-           sstime=float(nmind)
+           sstime=real(nmind,r_kind)
  
            tdiff=(sstime-gstime)*r60inv
 
diff --git a/src/gsi/read_ozone.f90 b/src/gsi/read_ozone.f90
index 43dd16d5c8..2ad95a858e 100644
--- a/src/gsi/read_ozone.f90
+++ b/src/gsi/read_ozone.f90
@@ -1056,7 +1056,7 @@ subroutine read_ozone(nread,ndata,nodata,jsatid,infile,gstime,lunout, &
            ozout(8,ndata)=usage1(k)          ! 
            ozout(9,ndata)=mlspres(k)         ! mls pressure in log(cb)
            ozout(10,ndata)=mlsozpc(k)        ! ozone mixing ratio precision in ppmv
-           ozout(11,ndata)=float(ipos(k))    ! pointer of obs level index in ozinfo.txt
+           ozout(11,ndata)=real(ipos(k),r_kind)    ! pointer of obs level index in ozinfo.txt
            ozout(12,ndata)=nloz              ! # of mls vertical levels
            ozout(nreal+1,ndata)=mlsoz(k)     ! ozone mixing ratio in ppmv
         end do
@@ -1220,7 +1220,7 @@ subroutine read_ozone(nread,ndata,nodata,jsatid,infile,gstime,lunout, &
          ozout(8,ndata)=usage1(k)          ! 
          ozout(9,ndata)=log(press(k))      ! ompslp pressure in log(cb)
          ozout(10,ndata)=omrstd(k)*ompslp_mult_fact   ! ozone mixing ratio precision in ppmv
-         ozout(11,ndata)=float(ipos(k))    ! pointer of obs level index in 
+         ozout(11,ndata)=real(ipos(k),r_kind)    ! pointer of obs level index in 
                                            ! ozinfo.txt
          ozout(12,ndata)=j !nloz              ! # of ompslp vertical levels
          ozout(13,ndata)=omr(k)            ! ozone mixing ratio in ppmv
diff --git a/src/gsi/read_pblh.f90 b/src/gsi/read_pblh.f90
index a7b7c066a2..1ff2cd2c23 100644
--- a/src/gsi/read_pblh.f90
+++ b/src/gsi/read_pblh.f90
@@ -343,7 +343,7 @@ subroutine read_pblh(nread,ndata,nodata,infile,obstype,lunout,twindin,&
 
 !       Add obs reference time, then subtract analysis time to get obs time relative to analysis
 
-        time_correction=float(minobs-minan)/60._r_kind
+        time_correction=real(minobs-minan,r_kind)/60._r_kind
 
       else
         time_correction=zero
diff --git a/src/gsi/read_prepbufr.f90 b/src/gsi/read_prepbufr.f90
index 304fa62590..eaece05451 100644
--- a/src/gsi/read_prepbufr.f90
+++ b/src/gsi/read_prepbufr.f90
@@ -1070,7 +1070,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
  
 !             Add obs reference time, then subtract analysis time to get obs time relative to analysis
  
-              time_correction=float(minobs-minan)*r60inv
+              time_correction=real(minobs-minan,r_kind)*r60inv
 
            else
               time_correction=zero
diff --git a/src/gsi/read_radar_wind_ascii.f90 b/src/gsi/read_radar_wind_ascii.f90
index 604d9d0eca..9d92699b6e 100644
--- a/src/gsi/read_radar_wind_ascii.f90
+++ b/src/gsi/read_radar_wind_ascii.f90
@@ -173,19 +173,19 @@ subroutine read_radar_wind_ascii(nread,ndata,nodata,infile,lunout,obstype,sis,hg
   end type radar
 
 !--Counters for diagnostics
- integer(i_kind) :: num_missing=0,numbadtime=0, &   !counts
-                    num_badtilt=0,num_badrange=0, &
-                    ibadazm=0 
+  integer(i_kind) :: num_missing=0,numbadtime=0, &   !counts
+                     num_badtilt=0,num_badrange=0, &
+                     ibadazm=0 
 
-integer(i_kind) :: ithin,zflag,nlevz,icntpnt,klon1,klat1,kk,klatp1,klonp1
-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
+  real(r_kind) :: rmesh,xmesh,zmesh,dx,dy,dx1,dy1,w00,w01,w10,w11
+  real(r_kind), allocatable, dimension(:) :: zl_thin
   real(r_kind),dimension(nsig):: hges,zges
   real(r_kind) sin2,termg,termr,termrg,zobs,height
   integer(i_kind) ntmp,iout,iiout,ntdrvr_thin2
   real(r_kind) crit1,timedif
   real(r_kind),parameter:: r16000 = 16000.0_r_kind
-logical :: luse
+  logical :: luse
   integer(i_kind) maxout,maxdata
   integer(i_kind),allocatable,dimension(:):: isort
 
@@ -201,7 +201,7 @@ subroutine read_radar_wind_ascii(nread,ndata,nodata,infile,lunout,obstype,sis,hg
   real(r_kind) :: thistiltr,selev0,celev0,thisrange,this_stahgt,thishgt                           
   real(r_kind) :: celev,selev,gamma,thisazimuthr,rlon0,t4dv, &
                   clat0,slat0,dlat,dlon,thiserr,thislon,thislat, &
-		  rlonloc,rlatloc,rlonglob,rlatglob,timeb,rad_per_meter
+                  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) :: rmins_an,rmins_ob                                                     
@@ -217,7 +217,7 @@ subroutine read_radar_wind_ascii(nread,ndata,nodata,infile,lunout,obstype,sis,hg
 
   type(radar),allocatable :: strct_in_vel(:,:)
 
-real(r_kind) :: mintilt,maxtilt,maxobrange,minobrange
+  real(r_kind) :: mintilt,maxtilt,maxobrange,minobrange
 
   integer(i_kind) :: thin_freq=1
 
@@ -228,10 +228,10 @@ subroutine read_radar_wind_ascii(nread,ndata,nodata,infile,lunout,obstype,sis,hg
 
   !-Check if radial velocity is in the convinfo file and extract necessary attributes 
   
- ithin=1 !number of obs to keep per grid box
- if(radar_no_thinning) then
-   ithin=-1
- endif
+  ithin=1 !number of obs to keep per grid box
+  if(radar_no_thinning) then
+    ithin=-1
+  endif
 
   errmax=-huge(errmax)
   errmin=huge(errmin)
@@ -241,13 +241,13 @@ subroutine read_radar_wind_ascii(nread,ndata,nodata,infile,lunout,obstype,sis,hg
      if(trim(obstype) == trim(ioctype(i)) .and. abs(icuse(i))== 1) then
         ikx=i 
         radartwindow=ctwind(ikx)*r60         !Time window units converted to minutes 
-	                                     !  (default setting for dbz within convinfo is 0.05 hours)
-	thiserr= 2_r_kind !1.75_r_kind !2_r_kind                     !Ob error (m/s) to use for radial velocity
-	exit                                 !Exit loop when finished with initial convinfo fields     
+                                             !  (default setting for dbz within convinfo is 0.05 hours)
+        thiserr= 2_r_kind !1.75_r_kind !2_r_kind                     !Ob error (m/s) to use for radial velocity
+        exit                                 !Exit loop when finished with initial convinfo fields     
      else if ( i==nconvtype ) then
         write(6,*) 'READ_RADAR_WIND_ASCII: ERROR - OBSERVATION TYPE IS NOT PRESENT IN CONVINFO OR USE FLAG IS ZERO'
-	write(6,*) 'READ_RADAR_WIND_ASCII: ABORTING read_radar_wind_ascii.f90 - NO VELOCITY OBS READ!'
-	return 	  
+        write(6,*) 'READ_RADAR_WIND_ASCII: ABORTING read_radar_wind_ascii.f90 - NO VELOCITY OBS READ!'
+        return   
      endif
   end do   
  
@@ -362,12 +362,12 @@ subroutine read_radar_wind_ascii(nread,ndata,nodata,infile,lunout,obstype,sis,hg
 
         obdate(1)=strct_in_vel(1,k)%year
         obdate(2)=strct_in_vel(1,k)%month  
-        obdate(3)=strct_in_vel(1,k)%day	 
+        obdate(3)=strct_in_vel(1,k)%day 
         obdate(4)=strct_in_vel(1,k)%hour   
         obdate(5)=strct_in_vel(1,k)%minute 
         call w3fs21(obdate,mins_ob)                             !mins_ob -integer number of mins snce 01/01/1978
-	rmins_ob=mins_ob                                        !convert to real number
-	rmins_ob=rmins_ob+(strct_in_vel(1,k)%second*r60inv)     !convert seconds to minutes and add to ob time
+        rmins_ob=mins_ob                                        !convert to real number
+        rmins_ob=rmins_ob+(strct_in_vel(1,k)%second*r60inv)     !convert seconds to minutes and add to ob time
  
       !-Comparison is done in units of minutes
       
@@ -377,70 +377,70 @@ subroutine read_radar_wind_ascii(nread,ndata,nodata,infile,lunout,obstype,sis,hg
         if(doradaroneob .and. (oneobradid /= strct_in_vel(1,k)%radid)) cycle tilts
 
         if(abs(timeb) > abs(radartwindow)) then
-	  numbadtime=numbadtime+1	  
-	  cycle tilts                           !If not in time window, cycle the loop
-	end if                  
+          numbadtime=numbadtime+1  
+          cycle tilts                           !If not in time window, cycle the loop
+        end if                  
       !--Time window check complete--!
 
         thistilt=strct_in_vel(1,k)%elev_angle
         if (thistilt <= maxtilt .and. thistilt >= mintilt) then 
      
           gates: do i=1,strct_in_vel(1,k)%num_gate,thin_freq  
-   	      thisrange=strct_in_vel(1,k)%fstgatdis + float(i-1)*strct_in_vel(1,k)%gateWidth
-	             
+              thisrange=strct_in_vel(1,k)%fstgatdis + real(i-1,r_kind)*strct_in_vel(1,k)%gateWidth
+             
              !-Check to make sure observations are within specified range 
 
-              if (thisrange <= maxobrange .and. thisrange >= minobrange) then	    
-	  	   
-	      azms: do j=1,strct_in_vel(1,k)%num_beam
-	   
-	           !-Check to see if this is a missing observation)
-		    nread=nread+1
-		    if ( strct_in_vel(1,k)%field(i,j) >= 999.0_r_kind ) then
-			  num_missing=num_missing+1
-	  	          cycle azms                        !No reason to process the ob if it is missing      	       
+              if (thisrange <= maxobrange .and. thisrange >= minobrange) then    
+     
+              azms: do j=1,strct_in_vel(1,k)%num_beam
+   
+                   !-Check to see if this is a missing observation)
+                    nread=nread+1
+                    if ( strct_in_vel(1,k)%field(i,j) >= 999.0_r_kind ) then
+                       num_missing=num_missing+1
+                       cycle azms                        !No reason to process the ob if it is missing      	       
                     end if
-		    		                      			
-		   !--Find observation height using method from read_l2bufr_mod.f90										       
-	         
-		    this_stahgt=strct_in_vel(1,k)%radhgt
+                        
+                   !--Find observation height using method from read_l2bufr_mod.f90										       
+         
+                    this_stahgt=strct_in_vel(1,k)%radhgt
                     aactual=rearth+this_stahgt                    
                     a43=four_thirds*aactual
                     thistiltr=thistilt*deg2rad
                     selev0=sin(thistiltr)
-                    celev0=cos(thistiltr)   		   
-		    b=thisrange*(thisrange+two*aactual*selev0)
+                    celev0=cos(thistiltr)      
+                    b=thisrange*(thisrange+two*aactual*selev0)
                     c=sqrt(aactual*aactual+b)
                     ha=b/(aactual+c)
                     epsh=(thisrange*thisrange-ha*ha)/(r8*aactual)
                     h=ha-epsh
-	            thishgt=this_stahgt+h 
+                    thishgt=this_stahgt+h 
                     height=thishgt
-		   !--Find observation location using method from read_l2bufr_mod.f90
-		 
-		   !-Get corrected tilt angle
-	            celev=celev0
-	            selev=selev0
-           	    celev=a43*celev0/(a43+h)
-		    selev=(thisrange*thisrange+h*h+two*a43*h)/(two*thisrange*(a43+h))
-	          
-		    gamma=half*thisrange*(celev0+celev)
-	         
+                   !--Find observation location using method from read_l2bufr_mod.f90
+
+                   !-Get corrected tilt angle
+                    celev=celev0
+                    selev=selev0
+                    celev=a43*celev0/(a43+h)
+                    selev=(thisrange*thisrange+h*h+two*a43*h)/(two*thisrange*(a43+h))
+          
+                    gamma=half*thisrange*(celev0+celev)
+         
                    !-Get earth lat lon of observation
-	         	 
+          
                     rlon0=deg2rad*strct_in_vel(1,k)%radlon
-	       	    clat0=cos(deg2rad*strct_in_vel(1,k)%radlat)
-		    slat0=sin(deg2rad*strct_in_vel(1,k)%radlat)		  
-		    thisazimuthr=(90.0_r_kind-strct_in_vel(1,k)%azim(j))*deg2rad   !Storing as 90-azm to
-		                                                                   ! be consistent with 
-										   ! read_l2bufr_mod.f90
-		    rad_per_meter=one/rearth
-		    rlonloc=rad_per_meter*gamma*cos(thisazimuthr)
+                    clat0=cos(deg2rad*strct_in_vel(1,k)%radlat)
+                    slat0=sin(deg2rad*strct_in_vel(1,k)%radlat)
+                    thisazimuthr=(90.0_r_kind-strct_in_vel(1,k)%azim(j))*deg2rad   !Storing as 90-azm to
+                                                                                   ! be consistent with 
+                                                                                   ! read_l2bufr_mod.f90
+                    rad_per_meter=one/rearth
+                    rlonloc=rad_per_meter*gamma*cos(thisazimuthr)
                     rlatloc=rad_per_meter*gamma*sin(thisazimuthr)
-		                  
-		    call invtllv(rlonloc,rlatloc,rlon0,clat0,slat0,rlonglob,rlatglob)
+                  
+                    call invtllv(rlonloc,rlatloc,rlon0,clat0,slat0,rlonglob,rlatglob)
                  
-		    thislat=rlatglob*rad2deg
+                    thislat=rlatglob*rad2deg
                     thislon=rlonglob*rad2deg 
 
                     if(doradaroneob) then
@@ -450,21 +450,21 @@ subroutine read_radar_wind_ascii(nread,ndata,nodata,infile,lunout,obstype,sis,hg
                      endif
 
 
-		    if(thislon>=r360) thislon=thislon-r360
+                    if(thislon>=r360) thislon=thislon-r360
                     if(thislon<zero ) thislon=thislon+r360
                  
-		   !-Convert back to radians                 
-         		       		       
-		    thislat = thislat*deg2rad
+                   !-Convert back to radians                 
+                
+                    thislat = thislat*deg2rad
                     thislon = thislon*deg2rad
-                 		 
-		    !find grid relative lat lon locations of earth lat lon
+                  
+                   !find grid relative lat lon locations of earth lat lon
                  
-		    call tll2xy(thislon,thislat,dlon,dlat,outside)
+                    call tll2xy(thislon,thislat,dlon,dlat,outside)
                     if (outside) cycle azms             !If observation is outside the domain
-		                                        ! then cycle, but don't increase range right away.
-							! Domain could be rectangular, so ob may be out of
-						        ! range at one end, but not the other.
+                                                        ! then cycle, but don't increase range right away.
+                                                        ! Domain could be rectangular, so ob may be out of
+                                                        ! range at one end, but not the other.
                 
                    if(regional .and. .not. fv3_regional) then
                        cosazm_earth=cos(thisazimuthr)
@@ -473,37 +473,37 @@ subroutine read_radar_wind_ascii(nread,ndata,nodata,infile,lunout,obstype,sis,hg
                        azm=atan2(sinazm,cosazm)
                     else
                        azm=thisazimuthr
-                    end if                                    		     					                   		   		   
-		   
-		   !--Do limited QC from read_radar.f90--!
-	            error = erradar_inflate*thiserr
+                    end if                                      
+
+                   !--Do limited QC from read_radar.f90--!
+                    error = erradar_inflate*thiserr
                     errmax=max(error,errmax)
                     if(thiserr>zero) errmin=min(error,errmin)                    
                     if(abs(azm)>r400) then
                        ibadazm=ibadazm+1
                        cycle azms
                     end if
-		    
-		    this_staid=strct_in_vel(1,k)%radid      !Via equivalence in declaration, value is propagated
-		    					    !  to rstation_id used below. 	    
-		    
-		    ! Get model terrain at radar station location
-		    ! If radar station is outside of grid, does not mean the 
-		    !    radar obs are outside the grid - therefore no need to
-		    !    cycle azms.
-		    
-		    radar_lon=deg2rad*strct_in_vel(1,k)%radlon
-		    radar_lat=deg2rad*strct_in_vel(1,k)%radlat
-		    call tll2xy(radar_lon,radar_lat,dlon_radar,dlat_radar,outside)
+
+                    this_staid=strct_in_vel(1,k)%radid      !Via equivalence in declaration, value is propagated
+                                                            !  to rstation_id used below. 	    
+    
+                    ! Get model terrain at radar station location
+                    ! If radar station is outside of grid, does not mean the 
+                    !    radar obs are outside the grid - therefore no need to
+                    !    cycle azms.
+    
+                    radar_lon=deg2rad*strct_in_vel(1,k)%radlon
+                    radar_lat=deg2rad*strct_in_vel(1,k)%radlat
+                    call tll2xy(radar_lon,radar_lat,dlon_radar,dlat_radar,outside)
                     call deter_zsfc_model(dlat_radar,dlon_radar,zsges)
-		    
-		    !  Determines land surface type based on surrounding land
+    
+                    !  Determines land surface type based on surrounding land
                     !    surface types
                      
-		    t4dv=timeb*r60inv
-		     
-	            call deter_sfc2(thislat,thislon,t4dv,idomsfc,skint,ff10,sfcr)
-		    
+                    t4dv=timeb*r60inv
+     
+                    call deter_sfc2(thislat,thislon,t4dv,idomsfc,skint,ff10,sfcr)
+    
 
 
 !####################       Data thinning       ###################
@@ -573,37 +573,37 @@ subroutine read_radar_wind_ascii(nread,ndata,nodata,infile,lunout,obstype,sis,hg
               isort(icntpnt)=iout
            endif
 
-		    cdata_all(1,iout) = error		                 ! wind obs error (m/s)
-		    cdata_all(2,iout) = dlon		                 ! grid relative longitude
-	            cdata_all(3,iout) = dlat 	                         ! grid relative latitude
-	     	    cdata_all(4,iout) = thishgt	                 ! obs absolute height (m)
-	     	    cdata_all(5,iout) = strct_in_vel(1,k)%field(i,j)    ! wind obs (m/s)
-	     	    cdata_all(6,iout) = azm                             ! azimuth angle (radians)
-	            cdata_all(7,iout) = t4dv       		         ! obs time (hour) - analysis relative
-	            cdata_all(8,iout) = ikx		                 ! type 	      
-	            cdata_all(9,iout) = thistiltr	                 ! tilt angle (radians)
-	            cdata_all(10,iout)= this_stahgt	                 ! station elevation (m)
-	            cdata_all(11,iout)= rstation_id	                 ! station id
-	            cdata_all(12,iout)= icuse(ikx)		         ! usage parameter
-	            cdata_all(13,iout)= idomsfc  	                 ! dominate surface type
-	            cdata_all(14,iout)= skint		                 ! skin temperature
-	            cdata_all(15,iout)= ff10		                 ! 10 meter wind factor
-	            cdata_all(16,iout)= sfcr		                 ! surface roughness
-	            cdata_all(17,iout)=thislon*rad2deg                  ! earth relative longitude (degrees)
-	            cdata_all(18,iout)=thislat*rad2deg                  ! earth relative latitude (degrees)
-	            cdata_all(19,iout)=thisrange/1000_r_kind		 ! range from radar in km (used to estimate beam spread)
-	            cdata_all(20,iout)=zsges		                 ! model elevation at radar site
-	            cdata_all(21,iout)=thiserr
-	            cdata_all(22,iout)=two                              ! Level 2 data
+                    cdata_all(1,iout) = error                            ! wind obs error (m/s)
+                    cdata_all(2,iout) = dlon                             ! grid relative longitude
+                    cdata_all(3,iout) = dlat                             ! grid relative latitude
+                    cdata_all(4,iout) = thishgt                          ! obs absolute height (m)
+                    cdata_all(5,iout) = strct_in_vel(1,k)%field(i,j)     ! wind obs (m/s)
+                    cdata_all(6,iout) = azm                              ! azimuth angle (radians)
+                    cdata_all(7,iout) = t4dv                             ! obs time (hour) - analysis relative
+                    cdata_all(8,iout) = ikx                              ! type 	      
+                    cdata_all(9,iout) = thistiltr                        ! tilt angle (radians)
+                    cdata_all(10,iout)= this_stahgt                      ! station elevation (m)
+                    cdata_all(11,iout)= rstation_id                      ! station id
+                    cdata_all(12,iout)= icuse(ikx)                       ! usage parameter
+                    cdata_all(13,iout)= idomsfc                          ! dominate surface type
+                    cdata_all(14,iout)= skint                            ! skin temperature
+                    cdata_all(15,iout)= ff10                             ! 10 meter wind factor
+                    cdata_all(16,iout)= sfcr                             ! surface roughness
+                    cdata_all(17,iout)=thislon*rad2deg                   ! earth relative longitude (degrees)
+                    cdata_all(18,iout)=thislat*rad2deg                   ! earth relative latitude (degrees)
+                    cdata_all(19,iout)=thisrange/1000._r_kind            ! range from radar in km (used to estimate beam spread)
+                    cdata_all(20,iout)=zsges                             ! model elevation at radar site
+                    cdata_all(21,iout)=thiserr
+                    cdata_all(22,iout)=two                               ! Level 2 data
 
                 if(doradaroneob .and. (cdata_all(5,iout) > -99_r_kind) ) exit volumes
 
                 end do azms  !j
               else
-	         num_badrange=num_badrange+1      !If outside acceptable range, increment
-	      end if   !Range check	
-		
-	   end do gates    !i
+                num_badrange=num_badrange+1      !If outside acceptable range, increment
+              end if   !Range check	
+
+           end do gates    !i
      
         else
            num_badtilt=num_badtilt+1           !If outside acceptable tilts, increment
diff --git a/src/gsi/read_saphir.f90 b/src/gsi/read_saphir.f90
index 3547f2cf5c..06e992b03d 100644
--- a/src/gsi/read_saphir.f90
+++ b/src/gsi/read_saphir.f90
@@ -360,10 +360,10 @@ subroutine read_saphir(mype,val_tovs,ithin,isfcalc,&
 
 
 ! compute look angle (panglr) and check against max angle
-!        panglr=(start+float(ifov-1)*step)*deg2rad
+!        panglr=(start+real(ifov-1,r_kind)*step)*deg2rad
 ! Use this calculation for now:
         step = .6660465
-        panglr = (42.96 - float(ifov-1)*step)*deg2rad
+        panglr = (42.96 - real(ifov-1,r_kind)*step)*deg2rad
 
         if(abs(lza)*rad2deg > MAX_SENSOR_ZENITH_ANGLE) then
           write(6,*)'READ_SAPHIR WARNING lza error ',lza,panglr
@@ -508,7 +508,7 @@ subroutine read_saphir(mype,val_tovs,ithin,isfcalc,&
 
 
 
-     crit1 = crit1 + rlndsea(isflg) + 10._r_kind*float(iskip) + 0.01_r_kind * abs(zz)
+     crit1 = crit1 + rlndsea(isflg) + 10._r_kind*real(iskip,r_kind) + 0.01_r_kind * abs(zz)
      call checkob(dist1,crit1,itx,iuse)
      if(.not. iuse)cycle ObsLoop
 
@@ -529,10 +529,10 @@ subroutine read_saphir(mype,val_tovs,ithin,isfcalc,&
      endif
 
 ! Re-calculate look angle
-!     panglr=(start+float(ifov-1)*step)*deg2rad
+!     panglr=(start+real(ifov-1,r_kind)*step)*deg2rad
 ! Use this calculation for now:
         step = .6660465
-        panglr = (42.96 - float(ifov-1)*step)*deg2rad
+        panglr = (42.96 - real(ifov-1,r_kind)*step)*deg2rad
 
 !     Load selected observation into data array
               
diff --git a/src/gsi/read_wcpbufr.f90 b/src/gsi/read_wcpbufr.f90
index 41dc36f3f6..f3daa5de43 100644
--- a/src/gsi/read_wcpbufr.f90
+++ b/src/gsi/read_wcpbufr.f90
@@ -430,7 +430,7 @@ subroutine read_wcpbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
  
 !             Add obs reference time, then subtract analysis time to get obs time relative to analysis
  
-              time_correction=float(minobs-minan)*r60inv
+              time_correction=real(minobs-minan,r_kind)*r60inv
 
            else
               time_correction=zero
diff --git a/src/gsi/reorg_metar_cloud.f90 b/src/gsi/reorg_metar_cloud.f90
index 2d947f2ef3..478d1ab363 100644
--- a/src/gsi/reorg_metar_cloud.f90
+++ b/src/gsi/reorg_metar_cloud.f90
@@ -264,9 +264,9 @@ subroutine reorg_metar_cloud(cdata,nreal,ndata,cdata_all,maxobs,ngrid)
              min_dist = 1.e10_r_kind
              do ic= 1,nsta_cld
                 ista = sta_cld(ic)
-                dist = (float(i1)-cdata(2,ista))*(float(i1)-cdata(2,ista))  &
-                      +(float(j1)-cdata(3,ista))*(float(j1)-cdata(3,ista))
-                if (dist < min_dist .and. dist < float(isprd2)) then
+                dist = (real(i1,r_kind)-cdata(2,ista))*(real(i1,r_kind)-cdata(2,ista))  &
+                      +(real(j1,r_kind)-cdata(3,ista))*(real(j1,r_kind)-cdata(3,ista))
+                if (dist < min_dist .and. dist < real(isprd2,r_kind)) then
                     min_dist = dist
                     ista_min = ista
                 end if
@@ -318,8 +318,8 @@ subroutine reorg_metar_cloud(cdata,nreal,ndata,cdata_all,maxobs,ngrid)
                       enddo
                       cdata_all(24,iout) = cdata_all(2,iout)   ! save observaion station i
                       cdata_all(25,iout) = cdata_all(3,iout)   ! save observaion station j
-                      cdata_all(2,iout) = float(i1)        ! grid index i
-                      cdata_all(3,iout) = float(j1)        ! grid index j
+                      cdata_all(2,iout) = real(i1,r_kind)        ! grid index i
+                      cdata_all(3,iout) = real(j1,r_kind)        ! grid index j
                       cdata_all(23,iout)= min_dist ! distance from station
                    endif
                 endif
diff --git a/src/gsi/rfdpar.f90 b/src/gsi/rfdpar.f90
index 79fa959bcb..b679ed6448 100644
--- a/src/gsi/rfdpar.f90
+++ b/src/gsi/rfdpar.f90
@@ -71,7 +71,7 @@ subroutine rfdpar1(be,rate,m)
   cof=zero
   cof(0)=one
   do i=1,m
-     cof(i)=half*cof(i-1)/float(i)
+     cof(i)=half*cof(i-1)/real(i,r_kind)
   enddo
 ! Locate the m roots of this polynomial:
   call zroots(cof,m,croot,polish)
diff --git a/src/gsi/satthin.F90 b/src/gsi/satthin.F90
index 2018d80be7..93f193014f 100644
--- a/src/gsi/satthin.F90
+++ b/src/gsi/satthin.F90
@@ -614,9 +614,9 @@ subroutine getsfc(mype,mype_io,use_sfc,use_sfc_any)
           jmax=nlat_sfc-2
           allocate(slatx(jmax),wlatx(jmax))
           call splat(idrt,jmax,slatx,wlatx)
-          dlon=two*pi/float(nlon_sfc)
+          dlon=two*pi/real(nlon_sfc,r_kind)
           do i=1,nlon_sfc
-             rlons_sfc(i)=float(i-1)*dlon
+             rlons_sfc(i)=real(i-1,r_kind)*dlon
           end do
           do i=1,(nlat_sfc-1)/2
              rlats_sfc(i+1)=-asin(slatx(i))
diff --git a/src/gsi/setupbend.f90 b/src/gsi/setupbend.f90
index 9bc856d67a..e82aa3dec9 100644
--- a/src/gsi/setupbend.f90
+++ b/src/gsi/setupbend.f90
@@ -331,9 +331,9 @@ subroutine setupbend(obsLL,odiagLL, &
 
 ! Intialize variables
   nsig_up=nsig+nsig_ext ! extend nsig_ext levels above interface level nsig
-  rsig=float(nsig)
+  rsig=real(nsig,r_kind)
   rdog=rd/grav
-  rsig_up=float(nsig_up)
+  rsig_up=real(nsig_up,r_kind)
   nobs_out=0
   hob_s_top=one
   mm1=mype+1
@@ -618,7 +618,7 @@ subroutine setupbend(obsLL,odiagLL, &
        ihob=hob
        k1=min(max(1,ihob),nsig)
        k2=max(1,min(ihob+1,nsig))
-       delz=hob-float(k1)
+       delz=hob-real(k1,r_kind)
        delz=max(zero,min(delz,one))
        trefges=tges_o(k1,i)*(one-delz)+tges_o(k2,i)*delz
        qrefges=qges_o(k1)*(one-delz)+qges_o(k2)*delz !Lidia
diff --git a/src/gsi/setupdw.f90 b/src/gsi/setupdw.f90
index 63c0df4b19..6ca68e4cae 100644
--- a/src/gsi/setupdw.f90
+++ b/src/gsi/setupdw.f90
@@ -298,7 +298,7 @@ subroutine setupdw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsa
   end if
 
   scale=one
-  rsig=float(nsig)
+  rsig=real(nsig,r_kind)
   mm1=mype+1
 
   call dtime_setup()
@@ -496,7 +496,7 @@ subroutine setupdw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsa
      dwwind=(ugesindw*sinazm+vgesindw*cosazm)*factw
 
      iz = max(1, min( int(dpres), nsig))
-     delz = max(zero, min(dpres - float(iz), one))
+     delz = max(zero, min(dpres - real(iz,r_kind), one))
 
      if (save_jacobian) then
         u_ind = getindex(svars3d, 'u')
diff --git a/src/gsi/setuplag.f90 b/src/gsi/setuplag.f90
index 0bad754a0c..01692d8164 100644
--- a/src/gsi/setuplag.f90
+++ b/src/gsi/setuplag.f90
@@ -171,7 +171,7 @@ subroutine setuplag(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diags
      allocate(cdiagbuf(nobs),rdiagbuf(nreal,nobs))
   end if
   scale=one
-  rsig=float(nsig)
+  rsig=real(nsig,r_kind)
   mm1=mype+1
 
   call dtime_setup()
diff --git a/src/gsi/setuplight.f90 b/src/gsi/setuplight.f90
index 040ef19bc6..b1118dd1f8 100644
--- a/src/gsi/setuplight.f90
+++ b/src/gsi/setuplight.f90
@@ -534,7 +534,7 @@ subroutine setuplight(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,light_di
 !            eps0     - guess value of lightning flash rate
 
      if(nobs_gbl > 0) then
-         eps=eps0*exp( (one/ float(nobs_gbl))*sum_gbl/(one+r0/w0) )
+         eps=eps0*exp( (one/ real(nobs_gbl,r_kind))*sum_gbl/(one+r0/w0) )
      else
          eps=eps0
      endif  !! if(nobs_gbl .gt. 0) then
diff --git a/src/gsi/setupoz.f90 b/src/gsi/setupoz.f90
index 7112e967ba..b16a33b414 100644
--- a/src/gsi/setupoz.f90
+++ b/src/gsi/setupoz.f90
@@ -1361,13 +1361,13 @@ subroutine setupozlev(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,&
 !    Check if observation above model top or below model surface
 
      rlow=max(sfcchk-dpres,zero)
-     rhgh=max(dpres-0.001_r_kind-float(nsig),zero)
+     rhgh=max(dpres-0.001_r_kind-real(nsig,r_kind),zero)
 
 !    calculate factor for error adjustment if too (high,low)
      ratio_errors=obserror/(obserror+1.0e6_r_kind*rhgh+four*rlow)
 
 !    Check to see if observations is above the top of the model
-     if (dpres > float(nsig)) then
+     if (dpres > real(nsig,r_kind)) then
          ratio_errors=zero
          obserror=1.0e6_r_kind
      end if
@@ -1379,7 +1379,7 @@ subroutine setupozlev(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,&
      call tintrp31(ges_oz,o3ges,dlat,dlon,dpres,dtime, &
        hrdifsig,mype,nfldsig)
      iz = max(1, min( int(dpres), nsig))
-     delz = max(zero, min(dpres - float(iz), one))
+     delz = max(zero, min(dpres - real(iz,r_kind), one))
      if (save_jacobian) then
         oz_ind = getindex(svars3d, 'oz')
         if (oz_ind < 0) then
diff --git a/src/gsi/setuppcp.f90 b/src/gsi/setuppcp.f90
index 8a6c8c0d80..7a89ccb6ca 100644
--- a/src/gsi/setuppcp.f90
+++ b/src/gsi/setuppcp.f90
@@ -416,8 +416,8 @@ subroutine setuppcp(obsLL,odiagLL,lunin,mype,aivals,nele,nobs,&
   elseif (amsu) then
      itype = 8
   endif
-  rterm1=one/float(nsig)
-  rterm2=one/float(nsig*(nsig-1))
+  rterm1=one/real(nsig,r_kind)
+  rterm2=one/real(nsig*(nsig-1),r_kind)
 
   call dtime_setup()
   do n = 1,nobs
diff --git a/src/gsi/setupq.f90 b/src/gsi/setupq.f90
index aa557b72c2..cebdeecd7b 100644
--- a/src/gsi/setupq.f90
+++ b/src/gsi/setupq.f90
@@ -674,7 +674,7 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
      error=one/(error*qsges)
 
      iz = max(1, min( int(dpres), nsig))
-     delz = max(zero, min(dpres - float(iz), one))
+     delz = max(zero, min(dpres - real(iz,r_kind), one))
 
 
      if (save_jacobian) then
diff --git a/src/gsi/setupref.f90 b/src/gsi/setupref.f90
index 752cce7e7b..4f240d756b 100644
--- a/src/gsi/setupref.f90
+++ b/src/gsi/setupref.f90
@@ -301,7 +301,7 @@ subroutine setupref(obsLL,odiagLL,lunin,mype,awork,nele,nobs,toss_gps_sub,is,ini
   ilate=15     ! index of earth relative latitude (degrees)
 
 ! Initialize variables
-  rsig=float(nsig)
+  rsig=real(nsig,r_kind)
   mm1=mype+1
 
 ! Check to see if required guess fields are available
@@ -967,11 +967,11 @@ subroutine setupref(obsLL,odiagLL,lunin,mype,awork,nele,nobs,toss_gps_sub,is,ini
               end do
            end if
 
-!          delz=dpres-float(k1)
+!          delz=dpres-real(k1,r_kind)
            kl=dpresl(i)
            k1l=min(max(1,kl),nsig)
            k2l=max(1,min(kl+1,nsig))
-           delz=dpresl(i)-float(k1l)
+           delz=dpresl(i)-real(k1l,r_kind)
            delz=max(zero,min(delz,one))
            my_head%jac_t(k1l)=my_head%jac_t(k1l)+termt(i)*(one-delz)
            my_head%jac_t(k2l)=my_head%jac_t(k2l)+termt(i)*delz
diff --git a/src/gsi/setupspd.f90 b/src/gsi/setupspd.f90
index 150799bf2c..2437ea63ce 100644
--- a/src/gsi/setupspd.f90
+++ b/src/gsi/setupspd.f90
@@ -504,7 +504,7 @@ subroutine setupspd(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diags
      spdges=sqrt(ugesin*ugesin+vgesin*vgesin)
 
      iz = max(1, min( int(dpres), nsig))
-     delz = max(zero, min(dpres - float(iz), one))
+     delz = max(zero, min(dpres - real(iz,r_kind), one))
 
      if (save_jacobian) then
 
diff --git a/src/gsi/setupt.f90 b/src/gsi/setupt.f90
index 5467a6dec9..d0ec421f06 100644
--- a/src/gsi/setupt.f90
+++ b/src/gsi/setupt.f90
@@ -515,7 +515,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
      if(netcdf_diag) call init_netcdf_diag_
   end if
   scale=one
-  rsig=float(nsig)
+  rsig=real(nsig,r_kind)
   mm1=mype+1
 
 !  rsli=isli
@@ -770,7 +770,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
                 hrdifsig,mype,nfldsig)
 
            iz = max(1, min( int(dpres), nsig))
-           delz = max(zero, min(dpres - float(iz), one))
+           delz = max(zero, min(dpres - real(iz,r_kind), one))
 
            if (save_jacobian) then
               t_ind = getindex(svars3d, 'tv')
@@ -792,7 +792,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
                 hrdifsig,mype,nfldsig)
 
            iz = max(1, min( int(dpres), nsig))
-           delz = max(zero, min(dpres - float(iz), one))
+           delz = max(zero, min(dpres - real(iz,r_kind), one))
 
            if (save_jacobian) then
               t_ind = getindex(svars3d, 'tsen')
diff --git a/src/gsi/setuptcp.f90 b/src/gsi/setuptcp.f90
index 3d13c5fe8e..68428175b1 100644
--- a/src/gsi/setuptcp.f90
+++ b/src/gsi/setuptcp.f90
@@ -255,9 +255,9 @@ subroutine setuptcp(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diags
         lb=dlat-tcp_box
         le=dlat+tcp_box
         do j=jb,je  
-           lj=float(j)
+           lj=real(j,r_kind)
            do l=lb,le
-              li=float(l)
+              li=real(l,r_kind)
               call tintrp2a11(ges_ps,psges,li,lj,dtime,hrdifsig,mype,nfldsig)
               if(pmin>psges)then
                  imin=l
diff --git a/src/gsi/setupw.f90 b/src/gsi/setupw.f90
index 784df1dfbe..6e653a9db0 100644
--- a/src/gsi/setupw.f90
+++ b/src/gsi/setupw.f90
@@ -677,7 +677,7 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
            hrdifsig,mype,nfldsig)
 
         iz = max(1, min( int(dpres), nsig))
-        delz = max(zero, min(dpres - float(iz), one))
+        delz = max(zero, min(dpres - real(iz,r_kind), one))
 
         if (save_jacobian) then
 
@@ -805,7 +805,7 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
            hrdifsig,mype,nfldsig)
 
         iz = max(1, min( int(dpres), nsig))
-        delz = max(zero, min(dpres - float(iz), one))
+        delz = max(zero, min(dpres - real(iz,r_kind), one))
 
         if (save_jacobian) then
            u_ind = getindex(svars3d, 'u')
diff --git a/src/gsi/setupwspd10m.f90 b/src/gsi/setupwspd10m.f90
index ad50c5b0c1..c702faaecf 100644
--- a/src/gsi/setupwspd10m.f90
+++ b/src/gsi/setupwspd10m.f90
@@ -449,7 +449,7 @@ subroutine setupwspd10m(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_d
            hrdifsig,mype,nfldsig)
  
         iz = max(1, min( int(dpres), nsig))
-        delz = max(zero, min(dpres - float(iz), one))
+        delz = max(zero, min(dpres - real(iz,r_kind), one))
 
         if (save_jacobian) then
 
diff --git a/src/gsi/sfcobsqc.f90 b/src/gsi/sfcobsqc.f90
index 8327b532e7..880e584b1b 100644
--- a/src/gsi/sfcobsqc.f90
+++ b/src/gsi/sfcobsqc.f90
@@ -1162,7 +1162,7 @@ subroutine get_wbinid(udbl,vdbl,nbins,ibin)
      endif
    else
      do n=1,nbins
-        if ( wdir >= float(n-1)*binwidth .and. wdir < float(n)*binwidth ) then 
+        if ( wdir >= real(n-1,r_kind)*binwidth .and. wdir < real(n,r_kind)*binwidth ) then 
             ibin=n
             exit
          endif
diff --git a/src/gsi/smoothzrf.f90 b/src/gsi/smoothzrf.f90
index 06e877ee68..dae656ca1b 100644
--- a/src/gsi/smoothzrf.f90
+++ b/src/gsi/smoothzrf.f90
@@ -81,7 +81,7 @@ subroutine frfhvo(p1,iv)
         do k=1,lat2
            l=int(rllat1(k,j))
            l2=min0(l+1,llmax)
-           dl2(k,j)=rllat1(k,j)-float(l)
+           dl2(k,j)=rllat1(k,j)-real(l,r_kind)
            dl1(k,j)=one-dl2(k,j)
         end do
      end do
diff --git a/src/gsi/ssmis_spatial_average_mod.f90 b/src/gsi/ssmis_spatial_average_mod.f90
index 24ac91d3f1..64dd5c6cf7 100644
--- a/src/gsi/ssmis_spatial_average_mod.f90
+++ b/src/gsi/ssmis_spatial_average_mod.f90
@@ -1551,7 +1551,7 @@ SUBROUTINE SFFTCB( X, N, M )
          END DO
          J = J + K
  104  CONTINUE
-      XT = 1.0 / FLOAT( N )
+      XT = 1.0 / real( N,r_kind )
       DO 99 I = 1, N
          X(I) = XT * X(I)
  99   CONTINUE
diff --git a/src/gsi/statsco.f90 b/src/gsi/statsco.f90
index ebef4a31e4..54f28ccfd9 100644
--- a/src/gsi/statsco.f90
+++ b/src/gsi/statsco.f90
@@ -119,7 +119,7 @@ subroutine statsco(stats_co,bwork,awork,ndata)
      if (iasim > 0) then
         svar = error_co(i)
         if (iuse_co(i)/=1) svar = -svar
-        rsum = one/float(iasim)
+        rsum = one/real(iasim,r_kind)
         icerr = nint(stats_co(2,i))
         do j=3,6   ! j=3=obs-mod(w_biascor)
                    ! j=4=(obs-mod(w_biascor))**2
@@ -145,7 +145,7 @@ subroutine statsco(stats_co,bwork,awork,ndata)
   do i=1,ndat
      if (idisplay(i)) then
         cpen=zero
-        if (icount_asim(i)>0) cpen=rpenal(i)/float(icount_asim(i))
+        if (icount_asim(i)>0) cpen=rpenal(i)/real(icount_asim(i),r_kind)
         write(iout_co,1115) jiter,dplat(i),dtype(i),ndata(i,2), &
              ndata(i,3),icount_asim(i),rpenal(i),cpen,qcpenal(i),iqccount_asim(i)
      endif
@@ -184,8 +184,8 @@ subroutine statsco(stats_co,bwork,awork,ndata)
         num(k)=nint(awork(5*nsig+k+100))
         rat=zero ; rat3=zero
         if(num(k) > 0) then
-           rat=awork(6*nsig+k+100)/float(num(k))
-           rat3=awork(3*nsig+k+100)/float(num(k))
+           rat=awork(6*nsig+k+100)/real(num(k),r_kind)
+           rat3=awork(3*nsig+k+100)/real(num(k),r_kind)
         end if
         ntot=ntot+num(k); o3plty=o3plty+awork(6*nsig+k+100)
         o3qcplty=o3qcplty+awork(3*nsig+k+100)
diff --git a/src/gsi/statsconv.f90 b/src/gsi/statsconv.f90
index 0da8606f24..3011fdefea 100644
--- a/src/gsi/statsconv.f90
+++ b/src/gsi/statsconv.f90
@@ -204,8 +204,8 @@ subroutine statsconv(mype,&
               rat1=zero
               rat2=zero
               if(num(k) > 0)then
-                 rat1=awork(4*nsig+k+100,i_uv)/float(num(k))
-                 rat2=awork(5*nsig+k+100,i_uv)/float(num(k))
+                 rat1=awork(4*nsig+k+100,i_uv)/real(num(k),r_kind)
+                 rat2=awork(5*nsig+k+100,i_uv)/real(num(k),r_kind)
               end if
               umplty=umplty+awork(4*nsig+k+100,i_uv)
               vmplty=vmplty+awork(5*nsig+k+100,i_uv)
@@ -218,8 +218,8 @@ subroutine statsconv(mype,&
               rat1=zero
               rat3=zero
               if(num(k) > 0)then
-                 rat1=(awork(4*nsig+k+100,i_uv)+awork(5*nsig+k+100,i_uv))/float(num(k))
-                 rat3=awork(3*nsig+k+100,i_uv)/float(num(k))
+                 rat1=(awork(4*nsig+k+100,i_uv)+awork(5*nsig+k+100,i_uv))/real(num(k),r_kind)
+                 rat3=awork(3*nsig+k+100,i_uv)/real(num(k),r_kind)
               end if
               uvqcplty=uvqcplty+awork(3*nsig+k+100,i_uv)
               write(iout_uv,240) 'w',num(k),k,awork(4*nsig+k+100,i_uv)+awork(5*nsig+k+100,i_uv), &
@@ -231,9 +231,9 @@ subroutine statsconv(mype,&
            write(iout_uv,925) 'wind',numgross,numfailqc
 !          Write statistics regarding penalties                   
            if(ntot > 0)then
-              tu=umplty/float(ntot)
-              tv=vmplty/float(ntot)
-              tuv=uvqcplty/float(ntot)
+              tu=umplty/real(ntot,r_kind)
+              tv=vmplty/real(ntot,r_kind)
+              tuv=uvqcplty/real(ntot,r_kind)
            end if
            if(numssm > 0)then
               tssm=awork(5,i_uv)/awork(6,i_uv)
@@ -286,8 +286,8 @@ subroutine statsconv(mype,&
               rat=zero
               rat3=zero
               if(num(k)>0) then
-                 rat=awork(6*nsig+k+100,i_gps)/float(num(k))
-                 rat3=awork(3*nsig+k+100,i_gps)/float(num(k))
+                 rat=awork(6*nsig+k+100,i_gps)/real(num(k),r_kind)
+                 rat3=awork(3*nsig+k+100,i_gps)/real(num(k),r_kind)
               end if
               ntot=ntot+num(k); gpsmplty=gpsmplty+awork(6*nsig+k+100,i_gps)
               gpsqcplty=gpsqcplty+awork(3*nsig+k+100,i_gps)
@@ -352,8 +352,8 @@ subroutine statsconv(mype,&
               rat=zero
               rat3=zero
               if(num(k) > 0)then
-                 rat=awork(5*nsig+k+100,i_q)/float(num(k))
-                 rat3=awork(3*nsig+k+100,i_q)/float(num(k))
+                 rat=awork(5*nsig+k+100,i_q)/real(num(k),r_kind)
+                 rat3=awork(3*nsig+k+100,i_q)/real(num(k),r_kind)
               end if
               qmplty=qmplty+awork(5*nsig+k+100,i_q)
               qqcplty=qqcplty+awork(3*nsig+k+100,i_q)
@@ -371,8 +371,8 @@ subroutine statsconv(mype,&
            numhgh      = nint(awork(3,i_q))
            write(iout_q,900) 'q',numhgh,numlow
            if(ntot > 0) then
-              tq=qmplty/float(ntot)
-              qctq=qqcplty/float(ntot)
+              tq=qmplty/real(ntot,r_kind)
+              qctq=qqcplty/real(ntot,r_kind)
            end if
         end if
 
@@ -414,8 +414,8 @@ subroutine statsconv(mype,&
            numfailqc=nint(awork(21,i_ps))
            write(iout_ps,925) 'psfc',numgross,numfailqc
            if(nump > 0)then
-              pw=awork(4,i_ps)/float(nump)
-              pw3=awork(22,i_ps)/float(nump)
+              pw=awork(4,i_ps)/real(nump,r_kind)
+              pw3=awork(22,i_ps)/real(nump,r_kind)
            end if
         end if
 
@@ -1116,8 +1116,8 @@ subroutine statsconv(mype,&
               num(k)=nint(awork(5*nsig+k+100,i_t))
               rat=zero ; rat3=zero
               if(num(k) > 0) then
-                 rat=awork(6*nsig+k+100,i_t)/float(num(k))
-                 rat3=awork(3*nsig+k+100,i_t)/float(num(k))
+                 rat=awork(6*nsig+k+100,i_t)/real(num(k),r_kind)
+                 rat3=awork(3*nsig+k+100,i_t)/real(num(k),r_kind)
               end if
               ntot=ntot+num(k); tmplty=tmplty+awork(6*nsig+k+100,i_t)
               tqcplty=tqcplty+awork(3*nsig+k+100,i_t)
@@ -1176,8 +1176,8 @@ subroutine statsconv(mype,&
               rat=zero
               rat3=zero
               if(num(k) > 0) then
-                 rat=awork(6*nsig+k+100,i_dw)/float(num(k))
-                 rat3=awork(3*nsig+k+100,i_dw)/float(num(k))
+                 rat=awork(6*nsig+k+100,i_dw)/real(num(k),r_kind)
+                 rat3=awork(3*nsig+k+100,i_dw)/real(num(k),r_kind)
               end if
               ntot=ntot+num(k)
               dwmplty=dwmplty+awork(6*nsig+k+100,i_dw)
@@ -1188,8 +1188,8 @@ subroutine statsconv(mype,&
            numgross=nint(awork(4,i_dw))
            numfailqc=nint(awork(21,i_dw))
            if(ntot > 0) then
-              tdw=dwmplty/float(ntot)
-              qctdw=dwqcplty/float(ntot)
+              tdw=dwmplty/real(ntot,r_kind)
+              qctdw=dwqcplty/real(ntot,r_kind)
            end if
            write(iout_dw,925) 'dw',numgross,numfailqc
            numlow       = nint(awork(2,i_dw))
@@ -1238,8 +1238,8 @@ subroutine statsconv(mype,&
               rat=zero
               rat3=zero
               if(num(k) > 0) then
-                 rat=awork(6*nsig+k+100,i_rw)/float(num(k))
-                 rat3=awork(3*nsig+k+100,i_rw)/float(num(k))
+                 rat=awork(6*nsig+k+100,i_rw)/real(num(k),r_kind)
+                 rat3=awork(3*nsig+k+100,i_rw)/real(num(k),r_kind)
               end if
               ntot=ntot+num(k)
               rwmplty=rwmplty+awork(6*nsig+k+100,i_rw)
@@ -1248,8 +1248,8 @@ subroutine statsconv(mype,&
                                               awork(3*nsig+k+100,i_rw),rat,rat3
            end do
            if(ntot > 0) then
-              trw=rwmplty/float(ntot)
-              qctrw=rwqcplty/float(ntot)
+              trw=rwmplty/real(ntot,r_kind)
+              qctrw=rwqcplty/real(ntot,r_kind)
            end if
            write(iout_rw,925) 'rw',numgross,numfailqc
            numlow       = nint(awork(2,i_rw))
@@ -1299,8 +1299,8 @@ subroutine statsconv(mype,&
               rat=zero
               rat3=zero
               if(num(k) > 0) then
-                 rat=awork(6*nsig+k+100,i_dbz)/float(num(k))
-                 rat3=awork(3*nsig+k+100,i_dbz)/float(num(k))
+                 rat=awork(6*nsig+k+100,i_dbz)/real(num(k),r_kind)
+                 rat3=awork(3*nsig+k+100,i_dbz)/real(num(k),r_kind)
               end if
               ntot=ntot+num(k)
               dbzmplty=dbzmplty+awork(6*nsig+k+100,i_dbz)
@@ -1309,8 +1309,8 @@ subroutine statsconv(mype,&
                                               awork(3*nsig+k+100,i_dbz),rat,rat3
            end do
            if(ntot > 0) then
-              tdbz=dbzmplty/float(ntot)
-              qctdbz=dbzqcplty/float(ntot)
+              tdbz=dbzmplty/real(ntot,r_kind)
+              qctdbz=dbzqcplty/real(ntot,r_kind)
            end if
            write(iout_dbz,925) 'dbz',numgross,numfailqc
            numlow       = nint(awork(2,i_dbz))
@@ -1360,8 +1360,8 @@ subroutine statsconv(mype,&
               rat=zero
               rat3=zero
               if(num(k) > 0) then
-                 rat=awork(6*nsig+k+100,i_fed)/float(num(k))
-                 rat3=awork(3*nsig+k+100,i_fed)/float(num(k))
+                 rat=awork(6*nsig+k+100,i_fed)/real(num(k),r_kind)
+                 rat3=awork(3*nsig+k+100,i_fed)/real(num(k),r_kind)
               end if
               ntot=ntot+num(k)
               fedmplty=fedmplty+awork(6*nsig+k+100,i_fed)
@@ -1370,8 +1370,8 @@ subroutine statsconv(mype,&
                                                awork(3*nsig+k+100,i_fed),rat,rat3
            end do
            if(ntot > 0) then
-              tfed=fedmplty/float(ntot)
-              qctfed=fedqcplty/float(ntot)
+              tfed=fedmplty/real(ntot,r_kind)
+              qctfed=fedqcplty/real(ntot,r_kind)
            end if
            write(iout_fed,925) 'fed',numgross,numfailqc
            numlow       = nint(awork(2,i_fed))
@@ -1419,8 +1419,8 @@ subroutine statsconv(mype,&
            write(iout_tcp,925) 'psfc',numgross,numfailqc
 
            if(nump > 0)then
-              pw=awork(4,i_tcp)/float(nump)
-              pw3=awork(22,i_tcp)/float(nump)
+              pw=awork(4,i_tcp)/real(nump,r_kind)
+              pw3=awork(22,i_tcp)/real(nump,r_kind)
            end if
         end if
 
@@ -1460,8 +1460,8 @@ subroutine statsconv(mype,&
               num(k)=nint(awork(6*nsig+k+100,i_lag))
               rat=zero ; rat3=zero
               if(num(k) > 0) then
-                 rat=awork(4*nsig+k+100,i_lag)/float(num(k))
-                 rat3=awork(3*nsig+k+100,i_lag)/float(num(k))
+                 rat=awork(4*nsig+k+100,i_lag)/real(num(k),r_kind)
+                 rat3=awork(3*nsig+k+100,i_lag)/real(num(k),r_kind)
               end if
               ntot=ntot+num(k); tmplty=tmplty+awork(4*nsig+k+100,i_lag)
               tqcplty=tqcplty+awork(3*nsig+k+100,i_lag)
diff --git a/src/gsi/statsoz.f90 b/src/gsi/statsoz.f90
index 069082d6b7..fb5f536914 100644
--- a/src/gsi/statsoz.f90
+++ b/src/gsi/statsoz.f90
@@ -101,7 +101,7 @@ subroutine statsoz(stats_oz,ndata)
      if (iasim > 0) then
         svar = error_oz(i)
         if (iuse_oz(i)/=1) svar = -svar
-        rsum = one/float(iasim)
+        rsum = one/real(iasim,r_kind)
         icerr = nint(stats_oz(2,i))
         do j=3,6   ! j=3=obs-mod(w_biascor)
                    ! j=4=(obs-mod(w_biascor))**2
@@ -127,7 +127,7 @@ subroutine statsoz(stats_oz,ndata)
   do i=1,ndat
      if (idisplay(i)) then
         cpen=zero
-        if (icount_asim(i)>0) cpen=rpenal(i)/float(icount_asim(i))
+        if (icount_asim(i)>0) cpen=rpenal(i)/real(icount_asim(i),r_kind)
         write(iout_oz,1115) jiter,dplat(i),dtype(i),ndata(i,2), &
              ndata(i,3),icount_asim(i),rpenal(i),cpen,qcpenal(i),iqccount_asim(i)
      endif
diff --git a/src/gsi/statspcp.f90 b/src/gsi/statspcp.f90
index e16f6dfbe5..ad79aed67d 100644
--- a/src/gsi/statspcp.f90
+++ b/src/gsi/statspcp.f90
@@ -210,7 +210,7 @@ subroutine statspcp(aivals,ndata)
      if (isum > 0 .and. display(is)) then
         rpen(is) = aivals(15,is)
         qcpen(is) = aivals(39,is)
-        rsum  = one/float(isum)
+        rsum  = one/real(isum,r_kind)
         icerr = nint(aivals(12,is))
         do j=13,16
            aivals(j,is)=aivals(j,is)*rsum
diff --git a/src/gsi/statsrad.f90 b/src/gsi/statsrad.f90
index 121761fa76..d42a53f7d6 100644
--- a/src/gsi/statsrad.f90
+++ b/src/gsi/statsrad.f90
@@ -120,7 +120,7 @@ subroutine statsrad(aivals,stats,ndata)
      if (iasim > 0) then
         svar = varch(i)
         if (iuse_rad(i) < 1) svar=-svar
-        rsum = one/float(iasim)
+        rsum = one/real(iasim,r_kind)
         icerr = nint(stats(2,i))
         do j=3,6   ! j=3=obs-mod(w_biascor)
                    ! j=4=(obs-mod(w_biascor))**2
diff --git a/src/gsi/stpjcmod.f90 b/src/gsi/stpjcmod.f90
index 2c811912a0..1cdabe60eb 100644
--- a/src/gsi/stpjcmod.f90
+++ b/src/gsi/stpjcmod.f90
@@ -871,7 +871,7 @@ subroutine stpjcpdry(rval,sval,pen,b,c,nbins)
   it=ntguessig
 
   dmass=zero_quad
-  rcon=one_quad/(two_quad*float(nlon))
+  rcon=one_quad/(two_quad*real(nlon,r_quad))
   mm1=mype+1
   return_now = .false.
   do n=1,nbins
diff --git a/src/gsi/support_2dvar.f90 b/src/gsi/support_2dvar.f90
index 35599ba548..0a83e96942 100644
--- a/src/gsi/support_2dvar.f90
+++ b/src/gsi/support_2dvar.f90
@@ -2465,19 +2465,19 @@ subroutine relocsfcob(rlon8,rlat8,cobtypein,cstationin,kxin)
      js=max(1,(jstart-jneighbour))
      je=min((jstart+jneighbour),ny)
 
-     ris=float(is)
-     rie=float(ie)
-     rjs=float(js)
-     rje=float(je)
+     ris=real(is,r_single)
+     rie=real(ie,r_single)
+     rjs=real(js,r_single)
+     rje=real(je,r_single)
 
      distmin=1.e+20_r_single
      lfound=.false.
 
      do j=1,npts
-      rj=rjs+float(j-1)*dy  
+      rj=rjs+real(j-1,r_single)*dy  
       if (rj > rje) cycle
       do i=1,npts
-         ri=ris+float(i-1)*dx  
+         ri=ris+real(i-1,r_single)*dx  
          if (ri > rie) cycle
 
          call bilinear_2d0(slmask,nx,ny,slmask0,rj,ri)
@@ -2655,7 +2655,7 @@ subroutine mkvalley_file
                  endif
               enddo
            enddo
-           hmean=hmean/max(1._r_single,float(ncount))
+           hmean=hmean/max(1._r_single,real(ncount,r_single))
 
            if ((hmax-hmin)>=hdiff0 .and. terrain(i,j)<hmean) & 
            valleys(i,j)=terrain(i,j)/hmean   
@@ -2744,8 +2744,8 @@ subroutine mkvalley_file
                     do ii=max(1,i-ijdel),min(nx,i+ijdel)
                        if (ii==i .and. jj==j) cycle
                        if (fldstd(ii,jj,1) < 1._r_single) then
-                          rminsq0=float((i-ii)*(i-ii)) + float((j-jj)*(j-jj))
-                          if (rminsq0 < rminsq .and. rminsq0 <= float(ijdel*ijdel)) then 
+                          rminsq0=real((i-ii)*(i-ii) + (j-jj)*(j-jj),r_single)
+                          if (rminsq0 < rminsq .and. rminsq0 <= real(ijdel*ijdel,r_single)) then 
                              rminsq=rminsq0
                              auxfld(i,j)=fldstd(ii,jj,1)  !Note that you are changing auxfld, not fldstd
                           endif
@@ -3515,8 +3515,8 @@ subroutine bilinear_2d0(rffcst,ix,jx,rfobs,xx,yy)
       i  = ifix(yy)
       j  = ifix(xx)
       
-      dx = xx - float(j)
-      dy = yy - float(i)
+      dx = xx - real(j,r_single)
+      dy = yy - real(i,r_single)
       dxm= 1.0_r_single-dx
       dym= 1.0_r_single-dy
  
diff --git a/src/gsi/tcv_mod.f90 b/src/gsi/tcv_mod.f90
index ab3e8c6372..90c886cd76 100644
--- a/src/gsi/tcv_mod.f90
+++ b/src/gsi/tcv_mod.f90
@@ -228,14 +228,14 @@ subroutine read_tcv_card(nums,storm,lucard,stswitch,slonfg,slatfg,centerid,stid,
           write (6,31) storm(i)
  
           if (storm(i)%tcv_lonew == 'W') then
-             slonfg(i) =  360._r_kind - float(storm(i)%tcv_lon)/10.0_r_kind
+             slonfg(i) =  360._r_kind - real(storm(i)%tcv_lon,r_kind)/10.0_r_kind
           else
-             slonfg(i) = float(storm(i)%tcv_lon)/10.0_r_kind
+             slonfg(i) = real(storm(i)%tcv_lon,r_kind)/10.0_r_kind
           endif
           if (storm(i)%tcv_latns == 'S') then
-             slatfg(i) = -one * float(storm(i)%tcv_lat)/10.0_r_kind
+             slatfg(i) = -one * real(storm(i)%tcv_lat,r_kind)/10.0_r_kind
           else
-             slatfg(i) = float(storm(i)%tcv_lat)/10.0_r_kind
+             slatfg(i) = real(storm(i)%tcv_lat,r_kind)/10.0_r_kind
           endif
         
           centerid(i) = storm(i)%tcv_center
diff --git a/src/gsi/tintrp2a.f90 b/src/gsi/tintrp2a.f90
index f675ae8b3f..d0a45124be 100644
--- a/src/gsi/tintrp2a.f90
+++ b/src/gsi/tintrp2a.f90
@@ -53,8 +53,8 @@ subroutine tintrp2a1(f,g,dx,dy,obstime,gridtime, &
   ix1=int(dx)
   iy1=int(dy)
   ix1=max(1,min(ix1,nlat))  
-  delx=dx-float(ix1)
-  dely=dy-float(iy1)
+  delx=dx-real(ix1,r_kind)
+  dely=dy-real(iy1,r_kind)
   delx=max(zero,min(delx,one))
   ix=ix1-istart(m1)+2
   iy=iy1-jstart(m1)+2
@@ -153,8 +153,8 @@ subroutine tintrp2a11(f,g,dx,dy,obstime,gridtime, &
   ix1=int(dx)
   iy1=int(dy)
   ix1=max(1,min(ix1,nlat))  
-  delx=dx-float(ix1)
-  dely=dy-float(iy1)
+  delx=dx-real(ix1,r_kind)
+  dely=dy-real(iy1,r_kind)
   delx=max(zero,min(delx,one))
   ix=ix1-istart(m1)+2
   iy=iy1-jstart(m1)+2
@@ -255,8 +255,8 @@ subroutine tintrp2a11_csln(f,g,gw,dx,dy,obstime,gridtime, &
   ix1=int(dx)
   iy1=int(dy)
   ix1=max(1,min(ix1,nlat))  
-  delx=dx-float(ix1)
-  dely=dy-float(iy1)
+  delx=dx-real(ix1,r_kind)
+  dely=dy-real(iy1,r_kind)
   delx=max(zero,min(delx,one))
   ix=ix1-istart(m1)+2
   iy=iy1-jstart(m1)+2
@@ -414,8 +414,8 @@ subroutine tintrp2a11_indx(dx,dy,obstime,gridtime, &
   ix1=int(dx)
   iy1=int(dy)
   ix1=max(1,min(ix1,nlat))
-  delx=dx-float(ix1)
-  dely=dy-float(iy1)
+  delx=dx-real(ix1,r_kind)
+  dely=dy-real(iy1,r_kind)
   delx=max(zero,min(delx,one))
   ix=ix1-istart(m1)+2
   iy=iy1-jstart(m1)+2
diff --git a/src/gsi/tintrp3.f90 b/src/gsi/tintrp3.f90
index c99830e649..39de029670 100644
--- a/src/gsi/tintrp3.f90
+++ b/src/gsi/tintrp3.f90
@@ -50,9 +50,9 @@ subroutine tintrp31(f,g,dx,dy,dz,obstime,gridtime,mype,nflds)
      iy1=int(dy)
      iz=int(dz)
      ix1=max(1,min(ix1,nlat)); iz=max(1,min(iz,nsig))  
-     delx=dx-float(ix1)
-     dely=dy-float(iy1)
-     delz=dz-float(iz)
+     delx=dx-real(ix1,r_kind)
+     dely=dy-real(iy1,r_kind)
+     delz=dz-real(iz,r_kind)
      delx=max(zero,min(delx,one)); delz=max(zero,min(delz,one))
      ix=ix1-istart(m1)+2
      iy=iy1-jstart(m1)+2
@@ -163,8 +163,8 @@ subroutine tintrp3(f,g,dx,dy,dz,obstime,gridtime,n,mype,nflds)
      ix1=int(dx)
      iy1=int(dy)
      ix1=max(1,min(ix1,nlat))
-     delx=dx-float(ix1)
-     dely=dy-float(iy1)
+     delx=dx-real(ix1,r_kind)
+     dely=dy-real(iy1,r_kind)
      delx=max(zero,min(delx,one))
      ix=ix1-istart(m1)+2
      iy=iy1-jstart(m1)+2
@@ -202,7 +202,7 @@ subroutine tintrp3(f,g,dx,dy,dz,obstime,gridtime,n,mype,nflds)
   do i=1,n
      iz=int(dz(i))
      iz=max(1,min(iz,nsig))  
-     delz=dz(i)-float(iz)
+     delz=dz(i)-real(iz,r_kind)
      delz=max(zero,min(delz,one))
      izp=min(iz+1,nsig)
      delzp=one-delz
diff --git a/src/gsi/wind_fft.f90 b/src/gsi/wind_fft.f90
index 449ceb120a..cc21843d8b 100644
--- a/src/gsi/wind_fft.f90
+++ b/src/gsi/wind_fft.f90
@@ -123,8 +123,8 @@ subroutine divvort_to_psichi(nx0,ny0,mmax0,nmax0,rld0,qg)
 
       pi = four*atan(one)
 
-      xmax = float(nx-1)
-      ymax = float(ny-1)
+      xmax = real(nx-1,r_kind)
+      ymax = real(ny-1,r_kind)
 !     write(6,*) ' in divvort_to_psichi: xmax,ymax,=',xmax,ymax
 !
 !==> compute trig tables for fft routines.
@@ -134,14 +134,14 @@ subroutine divvort_to_psichi(nx0,ny0,mmax0,nmax0,rld0,qg)
 !
 !==> set up wavenumbers used in fourier differentiation.
       do 100 i=1,nwavesx
-         rk(i) = two*pi*float(i-1)/xmax
+         rk(i) = two*pi*real(i-1,r_kind)/xmax
 100   continue
 
       ny2 = (ny/2)+1
       do 200 j=1,ny
          mm = j/ny2
          m = mm*ny+1
-         wavey = two*pi*float(j-m)/ymax
+         wavey = two*pi*real(j-m,r_kind)/ymax
          if (j<=nmax+1 .or. j>=ny-nmax+1) then
          indx = j-m+nmax+1
          indxy(indx) = j
diff --git a/src/gsi/write_fv3_spread.f90 b/src/gsi/write_fv3_spread.f90
index 4baef4c6bd..e1348c9449 100644
--- a/src/gsi/write_fv3_spread.f90
+++ b/src/gsi/write_fv3_spread.f90
@@ -227,10 +227,10 @@ subroutine write_fv3_enspread_ (grdin,filename,en_spread,ibin)
                          start = (/1/), count = (/grdin%nlon,grdin%nlat/)))
        ! levels
        do k=1,grdin%nsig
-         levsout(k) = float(k)
-         ilevsout(k) = float(k)
+         levsout(k) = real(k,r_kind)
+         ilevsout(k) = real(k,r_kind)
        end do
-       ilevsout(grdin%nsig+1) = float(grdin%nsig+1)
+       ilevsout(grdin%nsig+1) = real(grdin%nsig+1,r_kind)
        ! write to file
        call ncceck_enspread(nf90_put_var(ncid_out, levvarid, sngl(levsout), &
                          start = (/1/), count = (/grdin%nsig/)))
diff --git a/src/gsi/write_incr.f90 b/src/gsi/write_incr.f90
index 9adb150863..27a83d6e48 100644
--- a/src/gsi/write_incr.f90
+++ b/src/gsi/write_incr.f90
@@ -343,10 +343,10 @@ subroutine write_fv3_inc_ (grd,filename,mype_out,gfs_bundle,ibin)
                          start = (/1/), count = (/grd%nlon/)))
        ! levels
        do k=1,grd%nsig
-         levsout(k) = float(k)
-         ilevsout(k) = float(k)
+         levsout(k) = real(k,r_kind)
+         ilevsout(k) = real(k,r_kind)
        end do
-       ilevsout(grd%nsig+1) = float(grd%nsig+1)
+       ilevsout(grd%nsig+1) = real(grd%nsig+1,r_kind)
        ! write to file
        call nccheck_incr(nf90_put_var(ncid_out, levvarid, sngl(levsout), &
                          start = (/1/), count = (/grd%nsig/)))