c23456789012345678901234567890123456789012345678901234567890123456789012


      SUBROUTINE noising

c     ******************************************************************
c
c     This routine adds a random normal deviate to each component of the 
c     prior x_o 
c
c     ******************************************************************

c     -------------
c     Declarations
c     -------------

c     ''''''''''''''''''''''
      implicit none
      include 'commonvar.f'      
c     ''''''''''''''''''''''

c---- Local variables

      integer	i,j,k,idum

      integer	kc,kc1,kc2,kcr,kcomp

      real*8    n     (idim_n_o),
     >          n_rr  (idim_n_o),mean_n_r,var_n_r,
     >          cdf_sa(idim_n_o),
     >          cdf_no(idim_n_o),
     >          n_max

      real	gasdev

      logical	flagn(idim_n_o)


c---- NAG routine variables

      real*8	g01eaf

      character*1 tail

      integer	ifail


c---- External statement

c      external	g01eaf




c     ------------------------------------------------------------------
c     Adding noise to state variables (with zero mean and cte variance)
c     ------------------------------------------------------------------ 

c---- Initialization

      kcomp = 0


c---- Density 

      idum = -1

      do i=1,itot+1
      do j=1,jtot+1
      do k=1,ktot+1

         if (flag_r(i,j,k)) then

           kcomp    = kcomp + 1
c          ------------------------
           n(kcomp) = gasdev(idum)
c          ------------------------
           x(kcomp) = x(kcomp) + n(kcomp)*e_r(i,j,k)

         endif

      enddo
      enddo
      enddo


c---- Transport component U

      idum = -2

      do i=1,itot+1
      do j=1,jtot
      do k=1,ktot

         if (flag_U(i,j,k)) then

           kcomp    = kcomp + 1
c          ------------------------
           n(kcomp) = gasdev(idum)
c          ------------------------
           x(kcomp) = x(kcomp) + n(kcomp)*e_U(i,j,k)

        endif

      enddo
      enddo
      enddo

      
c---- Transport component V 

      idum = -3

      do i=1,itot
      do j=1,jtot+1
      do k=1,ktot

         if (flag_V(i,j,k)) then

           kcomp    = kcomp + 1
c          ------------------------
           n(kcomp) = gasdev(idum)
c          ------------------------
           x(kcomp) = x(kcomp) + n(kcomp)*e_V(i,j,k)

         endif

      enddo
      enddo
      enddo


c---- Transport component W 

      idum = -4

      do i=1  ,itot
      do j=1  ,jtot
      do k=1+1,ktot+1

         if (flag_W(i,j,k)) then

           kcomp    = kcomp + 1
c          ------------------------
           n(kcomp) = gasdev(idum)
c          ------------------------
           x(kcomp) = x(kcomp) + n(kcomp)*e_W(i,j,k)

         endif

      enddo
      enddo
      enddo


c---- Tracer

      idum = -5

      do i=1,itot+1
      do j=1,jtot+1
      do k=1,ktot+1

         if (flag_p(i,j,k)) then

           kcomp    = kcomp + 1
c          ------------------------
           n(kcomp) = gasdev(idum)
c          ------------------------
           x(kcomp) = x(kcomp) + n(kcomp)*e_p(i,j,k)

         endif

      enddo
      enddo
      enddo
 

c---- RETURN statement

      return




c     -----------------------------------------
c     CDF of the random normal deviate  
c
c     Note
c     ----
c     Comparison with the normal CDF permits
c     a subjective check
c     -----------------------------------------

c---- Ranking the deviates in ascending order

      do kc=1,kcomp
         flagn(kc) = .true.
      enddo

      n_max = -1.
      do kc=1,kcomp
         n_max = max(n_max,n(kc))
      enddo

      do kc1=1,kcomp

         n_rr(kc1) = n_max
         do kc2=1,kcomp
            if (n(kc2).lt.n_rr(kc1) .and. flagn(kc2)) then
              n_rr(kc1) = n(kc2)
              kcr       =   kc2
            endif
         enddo

         flagn(kcr) = .false.

      enddo


c---- CDF of deviates (no continuity correction)

      mean_n_r = 0.d0
      do kc=1,kcomp
         mean_n_r = mean_n_r + n_rr(kc)
      enddo
      mean_n_r = mean_n_r / dfloat(kcomp)
c--
      var_n_r = 0.d0
      do kc=1,kcomp
         var_n_r = var_n_r + (n_rr(kc)-mean_n_r)**2
      enddo
      var_n_r = var_n_r / dfloat(kcomp-1) 
c--
      do kc=1,kcomp
         n_rr(kc) = (n_rr(kc)-mean_n_r) / sqrt(var_n_r)
      enddo
c--
      do kc1=1,kcomp
         cdf_sa(kc1) = 0.d0
         do kc2=1,kcomp
            if (n_rr(kc2).le.n_rr(kc1)) cdf_sa(kc1) = cdf_sa(kc1) + 1.d0  
         enddo
         cdf_sa(kc1) = cdf_sa(kc1) / dfloat(kcomp)
      enddo


c---- CDF of the standard normal distribution 

      ifail = 0

      do kc=1,kcomp

c        ---------------------------------------
c         cdf_no(kc) = g01eaf('L',n_rr(kc),ifail)
c        ---------------------------------------

         if (ifail.ne.0) then
           write(*,*)
           write(*,*) 'ifail=',ifail,' in noising.f'
           write(*,*)
           stop 'STOP'
         endif

      enddo


c---- Writing 

      open (19,file='test')
      do kc=1,kcomp
         write(19,*) n_rr(kc),cdf_sa(kc),
     >			     cdf_no(kc)
      enddo
      close (19)




c     ----------------------
c     End of the subroutine
c     ----------------------

      return

      end

c     ******************************************************************


      REAL FUNCTION gasdev(idum)

c     ******************************************************************
c
c     This function returns a normally ditributed deviate with zero mean
c     and unit variance, using ran1(idum) as the source of uniform deviate
c
c
c     Method: Press et al., Numerical Recipes in Fortran, 1992
c     ------
c    
c     ******************************************************************

c     -------------
c     Declarations
c     -------------

c     ''''''''''''''
      implicit none
c     ''''''''''''''

c---- Argument variable

      integer	idum


c---- Local variables

      integer   iset

      real      fac,gset,rsq,v1,v2,ran1

      save      iset,gset

      data      iset/0/


c     ----------------------------------------------
c     Generation of the sequence of normal deviates
c     ----------------------------------------------

c++
      if (iset.eq.0) then
c++
1       v1  = 2.*ran1(idum)-1.
        v2  = 2.*ran1(idum)-1.
        rsq = v1**2+v2**2
        if (rsq.ge.1..or.rsq.eq.0.) goto 1
        fac    = sqrt(-2.*log(rsq)/rsq)
        gset   = v1*fac
        gasdev = v2*fac
        iset   = 1
c++
      else
c++
        gasdev = gset
        iset   = 0
c++
      endif
c++


c     --------------------
c     End of the function
c     --------------------

      return

      end

c     ******************************************************************


      REAL FUNCTION ran1(idum)

c     ******************************************************************
c
c     This function returns a uniform random deviate between 0.0 and 1.0
c     Set idum to any negative integer to initialize. Do not alter idum
c     between successive deviates in a sequence.
c
c     Method: Press et al., Numerical Recipes in FORTRAN, 1992
c     ------
c     
c     ******************************************************************

c---- Declarations

c     ''''''''''''''     
      implicit none
c     ''''''''''''''


c---- Argument variable

      integer	idum


c---- Local variables

      integer   ia,im,iq,ir,ntab,ndiv

      real      am,eps,rnmx

      parameter (ia=16807,im=2147483647,iq=127773,ir=2836)

      parameter (ntab=32,ndiv=1+(im-1)/ntab)

      parameter (am=1./im,eps=1.2e-7,rnmx=1.-eps)

      integer   j,k,iv(ntab),iy

      save      iv,iy

      data iv /ntab*0/, iy /0/


c---- Generation of a sequence of random numbers

c++
      if (idum.le.0 .or. iy.eq.0) then
c++
        idum = max(-idum,1)
        do j=ntab+8,1,-1
           k    = idum / iq
           idum = ia*(idum-k*iq) - ir*k
           if (idum.lt.   0) idum  = idum + im
           if (j   .le.ntab) iv(j) = idum
        enddo
        iy = iv(1)
c++
      endif
c++

      k    = idum / iq
      idum = ia*(idum-k*iq) - ir*k
      if (idum.lt.0) idum = idum + im
      j     = 1 + iy/ndiv
      iy    = iv(j)
      iv(j) = idum
      ran1  = min(am*iy,rnmx)



c---- End of the function

      return

      end

c     ***********************************************************************



