c=========================================================================
c
c NAME: grid_moment_cs2
c
c PURPOSE:  Computes mean and standard deviation from summed arrays.
c
c INPUTS: name   data type   dimension   description
c         ---------------------------------------------------------
c         nlon   integer     scalar      number of longitude points
c         nlat   integer     scalar      number of latitude points
c         sum    real        3D          summed value at given gridpoint
c         sumsq  real        3D          sum**2 value at given gridpoint.(float)
c         count  real        3D          total number of points.(Int*4)
c
c OUTPUTS:
c         imn    integer*2   3D          mean value
c         isd    integer*2   3D          standard deviation
c         ict    integer*2   3D          total number of points
c
c         3D = nlon x nlat x nparm
c
c EXTERNALS:
c            hirs_grid.inc
c            iscale.f
c
c AUTHOR:
c         Darren Jackson CIRES/ETL April 2002
c           - version 2, January 2004, takes into account asc/des passes 
c                        for clear-sky grid data
c
c============================================================================
      subroutine grid_moment_cs2(nlon,nlat,sum,sumsq,count,
     .     imn,isd,ict)

      implicit none
      include 'hirs_grid.inc'
      
      integer nlon,nlat,i,j,k

      integer*2 ict(nlon,nlat,node),imn(nlon,nlat,node)
      integer*2 isd(nlon,nlat,node),iscale

      integer count(nlon,nlat,node)
      real sum(nlon,nlat,node),sumsq(nlon,nlat,node),mn,sd,tmp1,tmp2
      real offset_mn,sfactor_mn,offset_sd,sfactor_sd

      do k=1,node
         do j=1,nlat
            do i=1,nlon
               if(count(i,j,k) .gt. 0) then
                  mn=sum(i,j,k)/count(i,j,k)
                  tmp1=sumsq(i,j,k)/count(i,j,k)
                  tmp2=(sum(i,j,k)/count(i,j,k))**2.
                  if(tmp1 .ge. tmp2) then
                     sd=sqrt(tmp1-tmp2)
                  else
                     sd=0.
                  endif
               else
                  mn=imiss
                  sd=imiss
               endif
            
               offset_mn = 100.
               sfactor_mn = 100.
               offset_sd = 0.
               sfactor_sd = 100.

               imn(i,j,k)=iscale(mn,offset_mn,sfactor_mn)
               isd(i,j,k)=iscale(sd,offset_sd,sfactor_sd)
               ict(i,j,k)=count(i,j,k)

            enddo
         enddo
      enddo

      return
      end