Merge remote-tracking branch 'emc/develop' into CADS_for_Andrew

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
 

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

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

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