program random
! **********************************************************************
! *              .      .    .
! * PROGRAM: random
! *
! * PRGMMR:     M. ZUPANSKI, D. ZUPANSKI     ORG: CIRA/CSU   DATE: 2003-08-14
! *
! * ABSTRACT:   THIS SUBROUTINE CONVERTS CONTROL (FCST) VARIABLES FROM
! *             1-DIMENSIONAL TO N-DIMENSIONAL VARIABLE, 
! *             AND WRITES OUT THE N-DIMENSIONAL VARIABLE.
! *
! * PROGRAM LOG:
! *
! * 08/14/2003    ..... M. ZUPANSKI, D. ZUPANSKI
! * 09/12/2003    ..... D. ZUPANSKI
! * 09/18/2004    ..... M. ZUPANSKI: Add KPZ equation to 'smooth.F'
! *
! **********************************************************************

#ifdef MPI_USE
  include "mpif.h"
#endif

      integer,parameter :: icntrl=20           ! input control variable
      integer,parameter::idev=21              ! input standard deviation
      integer,parameter :: i_ic=51             ! initial conditions file
      integer,parameter :: i_param=52          ! empirical parameter file
      integer,parameter :: i_bias=53           ! starting bias file
      integer,parameter :: iptrb=121           ! output perturbed (ensemble) cntrl vrbl

      INCLUDE 'define_var_list_type.h' 
      INCLUDE 'define_cvar_type.h' 

      TYPE (var_list_type), dimension(:), allocatable :: cvar_list
      TYPE (cvar_type), dimension(:), allocatable :: cvar
  
!--- random numbers
       integer,parameter::imax=1e6
       integer :: k
       integer,dimension(:),allocatable :: seed
       real,dimension(:),allocatable :: random_x
       real,dimension(:),allocatable :: stddev
       real,dimension(:),allocatable :: stddev1d
!      real,dimension(:,:),allocatable :: stddev2d
       real,dimension(:,:,:),allocatable :: stddev3d
!      real,dimension(:,:,:,:),allocatable :: stddev4d
!-----
   character*8 fileout
           
       integer,allocatable :: jlen(:),jdisp(:)
       integer :: N_LOC,jsta,jend
       integer :: NPROC,MPIRANK
       integer :: NENS
       integer :: NENS_START
       integer :: NALL

      real, dimension(:), allocatable :: x,x_cntrl
      real, dimension(:), allocatable :: cvar1d
      real, dimension(:,:), allocatable :: cvar2d
      real, dimension(:,:,:), allocatable :: cvar3d
      real, dimension(:,:,:,:), allocatable :: cvar4d
      real                                  :: cvarparm
      real                                  :: correl
      integer, dimension(1) :: NNXP,NNYP,NNZP,NNSOIL,NNSNOW
      integer, dimension(:), allocatable :: unitbias
      integer               :: cvar_max
      integer               :: ubias,unitbfile,nbias,num_bias
      character*20          :: filename
      character*20          :: file_ic,file_param,file_bias 
      character*9           :: cv_name
      logical               :: open_ic,open_param,open_bias

!------------------------------
      NAMELIST /MODEL_DIMENSION/ NNXP,NNYP,NNZP,NNSOIL,NNSNOW
      NAMELIST /ENSEMBLE_SIZE/ NENS_START,NENS
!------------------------------

#ifdef MPI_USE
!
!     START MPI
!
      CALL MPI_INIT(IERR)
      IF (IERR .NE. MPI_SUCCESS) STOP 'FAILED TO INIT MPI'
      CALL MPI_COMM_SIZE(MPI_COMM_WORLD, NPROC, IERR)
      CALL MPI_COMM_RANK(MPI_COMM_WORLD, MPIRANK, IERR)
      CALL MPI_BARRIER(MPI_COMM_WORLD,IERR)
#else
      write(*,*) "This is a NO MPI run"
      MPIRANK=0
      NPROC=1
#endif

       if(MPIRANK.eq.0) then
       write(*,*) "start random"
       end if

!--- get model dimensions 
      write(filename,200)
 200  format('model_dimension.name')
      CLOSE(500)
      OPEN(UNIT=500,FILE=filename,FORM='FORMATTED',IOSTAT=IER)
      IF(IER.NE.0) WRITE(3,*)' 500 OPEN UNIT ERROR IER=',IER
      READ(500,MODEL_DIMENSION)  
      CLOSE(500,STATUS='KEEP')
       if(MPIRANK.eq.0) then
       write(*,*) "NNXP=",NNXP(1)," NNYP=",NNYP(1)," NNZP=",NNZP(1)
       end if

      INCLUDE 'cntrl_vrbl_list.h' 

             if(MPIRANK.eq.0) then
      write(*,*) "max_num_of_cntrl_vrbls=",max_num_of_cntrl_vrbls
       do ncv=1,max_num_of_cntrl_vrbls
      write(*,*) "ndim ",cvar_list(ncv)%ndim
      write(*,*) "start1 ",cvar_list(ncv)%start_index(1)
      write(*,*) "start2 ",cvar_list(ncv)%start_index(2)
      write(*,*) "start3 ",cvar_list(ncv)%start_index(3)
      write(*,*) "start4 ",cvar_list(ncv)%start_index(4)
      write(*,*) "end1 ",cvar_list(ncv)%end_index(1)
      write(*,*) "end2 ",cvar_list(ncv)%end_index(2)
      write(*,*) "end3 ",cvar_list(ncv)%end_index(3)
      write(*,*) "end4 ",cvar_list(ncv)%end_index(4)
      write(*,*) "name ",cvar_list(ncv)%name
      write(*,*) "description ",cvar_list(ncv)%description
       end do
             end if

!-- read ensemble size
                                                                                               
       REWIND 15
       READ(15,ENSEMBLE_SIZE)
                                                                                               
          if(MPIRANK.eq.0) then
          write(*,*) "random: NENS_START, NENS=",NENS_START,NENS
          end if
 
!-- Total number of degrees of freedom
     NALL=NENS-NENS_START+1


      CLOSE(103)
      OPEN(103,FILE='cntrl_vrbl',FORM='FORMATTED',IOSTAT=IER)
      IF(IER.NE.0)WRITE(*,*)' 103 OPEN UNIT ERROR IER=',IER
      READ(103,*) cvar_max
      allocate(cvar(1:cvar_max))
      open_ic=.FALSE.
      open_param=.FALSE.
      open_bias=.FALSE.
       do i=1,cvar_max
      READ(103,*) cvar(i)
          if(MPIRANK.eq.0) then
           write(*,*) cvar(i)
            end if
         if(cvar(i)%ic) then
         open_ic=.TRUE.
         end if
           if(cvar(i)%param) then
           open_param=.TRUE.
           end if
              if(cvar(i)%bias) then
              open_bias=.TRUE.
              num_bias=cvar(i)%num_bias
              end if
!     READ(103,*) name,ic,param,bias,num_bias !! read in the model
!     READ(103,*) cvar(i)%name,cvar(i)%ic,cvar(i)%param,cvar(i)%bias,cvar(i)%num_bias
!     write(*,*) cvar(i)%name,cvar(i)%ic,cvar(i)%param,cvar(i)%bias,cvar(i)%num_bias
       end do
      CLOSE(103,STATUS='KEEP')

!-- read 1-d control variable
 
      CLOSE(icntrl)
      OPEN(icntrl,file='x1d_cntrl',status='unknown',form='unformatted')
      READ(icntrl) idim_1d
        if(MPIRANK.eq.0) then
        write(*,*) "x dimension=",idim_1d
        end if
      allocate(x_cntrl(1:idim_1d))
      READ(icntrl) x_cntrl
      CLOSE(icntrl,status='KEEP')


!-- read standard deviation
                                                                               
      CLOSE(idev)
      OPEN(idev,file='std_dev',status='unknown',form='unformatted')
      READ(idev) jdim
        if(MPIRANK.eq.0) then
        write(*,*) "stdev dimension=",jdim
        end if
      allocate(stddev(1:idim_1d))
      READ(idev) stddev
      CLOSE(idev,status='KEEP')

!====================================================================
!==
!==  define local dimensions and indexes
!==
       allocate(jdisp(0:NPROC-1)) ; allocate(jlen(0:NPROC-1))
                                                                                               
            do irank=0,NPROC-1
      CALL para_range(NENS_START,NENS,NPROC,irank,jsta,jend)
            jlen(irank)=jend-jsta+1
            jdisp(irank)=jsta-NENS_START
            end do
#ifdef MPI_USE
      CALL MPI_BARRIER(MPI_COMM_WORLD,IERR)
#endif
      CALL para_range(NENS_START,NENS,NPROC,MPIRANK,jsta,jend)
                                                                                               
!==
         N_LOC=jlen(MPIRANK)
!==
!
         if(MPIRANK.eq.0) then
           do nn=0,NPROC-1
           write(*,*) "nn=",nn," jdisp=",jdisp(nn)," jlen=",jlen(nn)
           end do
          end if
!=================================================================
!--  SETUP INITIAL SEED -------------------------------
       
       allocate(random_x(1:idim_1d))
        do ii=1,1000
         call random_number(random_x(1:idim_1d))
          end do
       deallocate(random_x)
       
#ifdef MPI_USE
      CALL MPI_BARRIER(MPI_COMM_WORLD,IERR)
#endif
       
       call random_seed(size=k)
        allocate(seed(1:k))
         call random_seed(get=seed(1:k))
       
!-- change seed for each CPU --------------------------
      seed(1:k)=seed(1:k)+mpirank*10
       
!-- save new seed --------------------------------------
      call random_seed(put=seed(1:k))
       
      deallocate(seed)
!===================================================
      allocate(x(1:idim_1d))

      DO N=jsta,jend                  !!! loop over perturbations (ensembles)

!-- Generate default perturbations as Gauss random noise----------
!-- Gauss random number generator ------------
       allocate(random_x(1:idim_1d))
      mode2=MOD(idim_1d,2)
      Nsave=idim_1d/2 + mode2
      do l=1,Nsave
      call Box_Muller_polar(y1,y2)
!----------------------------------------------
      ii_tot=0
          do while (abs(y1).gt.3..or.abs(y2).gt.3.)
      ii_tot=ii_tot+1
      if(ii_tot.gt.10) exit
      call Box_Muller_polar(yy1,yy2)
          if(abs(yy1).le.3.) then
          y1=yy1
          end if
            if(abs(yy2).le.3.) then
            y2=yy2
            end if
         end do
                                                                               
      random_x(l)=y1
      random_x(Nsave+l-mode2)=y2
      end do
!--------------------------------------------
     x(:)=x_cntrl(:) + random_x(:)*stddev(:)
!--------------------------------------------
       deallocate(random_x)

        icount=0
        icount_out=0
        nvar=0
       do i=1,cvar_max
        do j=1,max_num_of_cntrl_vrbls
        if(cvar(i)%name.eq.cvar_list(j)%name) then
        nvar=nvar+1
        correl=cvar_list(j)%correl

         i1s=cvar_list(j)%start_index(1)
         i2s=cvar_list(j)%start_index(2)
         i3s=cvar_list(j)%start_index(3)
         i4s=cvar_list(j)%start_index(4)
         i1e=cvar_list(j)%end_index(1)  
         i2e=cvar_list(j)%end_index(2)  
         i3e=cvar_list(j)%end_index(3)  
         i4e=cvar_list(j)%end_index(4)  

!-- initial conditions
         if(cvar(i)%ic) then
         cv_name=cvar(i)%name
          if(cvar_list(j)%ndim.eq.4) then
          allocate(cvar4d(i1s:i1e,i2s:i2e,i3s:i3e,i4s:i4e))
!         allocate(stddev4d(i1s:i1e,i2s:i2e,i3s:i3e,i4s:i4e))
            do i4=i4s,i4e
             do i3=i3s,i3e
              do i2=i2s,i2e
               do i1=i1s,i1e
                icount=icount+1
                cvar4d(i1,i2,i3,i4)=x_cntrl(icount)
!               stddev4d(i1,i2,i3,i4)=stddev(icount)
               end do 
              end do 
             end do 
            end do 

!     call perturb_4d(cv_name,i1s,i1e,i2s,i2e,i3s,i3e,i4s,i4e, &
!                 stddev4d,correl,cvar4d,MPIRANK)
!         deallocate(stddev4d)
! For the moment we use default perturbations

            do i4=i4s,i4e
             do i3=i3s,i3e
              do i2=i2s,i2e
               do i1=i1s,i1e
                icount_out=icount_out+1
!               x(icount_out)=cvar4d(i1,i2,i3,i4)
               end do 
              end do 
             end do 
            end do 

          deallocate(cvar4d)
           else if(cvar_list(j)%ndim.eq.3) then
           allocate(cvar3d(i1s:i1e,i2s:i2e,i3s:i3e))
           allocate(stddev3d(i1s:i1e,i2s:i2e,i3s:i3e))
             do i3=i3s,i3e
              do i2=i2s,i2e
               do i1=i1s,i1e
                icount=icount+1
                cvar3d(i1,i2,i3)=x_cntrl(icount)
                stddev3d(i1,i2,i3)=stddev(icount)
               end do 
              end do 
             end do 

    write(*,*) " min max stddev3d=",minval(stddev3d),maxval(stddev3d)

   call perturb_3d (cv_name,i1s,i1e,i2s,i2e,i3s,i3e,  &
                    stddev3d,correl,cvar3d,MPIRANK)

          deallocate(stddev3d)

             do i3=i3s,i3e
              do i2=i2s,i2e
               do i1=i1s,i1e
                icount_out=icount_out+1
                x(icount_out)=cvar3d(i1,i2,i3)
               end do 
              end do 
             end do 

           deallocate(cvar3d)
            else if(cvar_list(j)%ndim.eq.2) then
            allocate(cvar2d(i1s:i1e,i2s:i2e))
!           allocate(stddev2d(i1s:i1e,i2s:i2e))
              do i2=i2s,i2e
               do i1=i1s,i1e
                icount=icount+1
                cvar2d(i1,i2)=x_cntrl(icount)
!               stddev2d(i1,i2,i3)=stddev(icount)
               end do 
              end do 

!  call perturb_2d (cv_name,i1s,i1e,i2s,i2e,stddev2d,correl,cvar2d,MPIRANK)

!         deallocate(stddev2d)

              do i2=i2s,i2e
               do i1=i1s,i1e
                icount_out=icount_out+1
!               x(icount_out)=cvar2d(i1,i2)
               end do 
              end do 

            deallocate(cvar2d)
             else if(cvar_list(j)%ndim.eq.1) then
             allocate(cvar1d(i1s:i1e))
             allocate(stddev1d(i1s:i1e))
               do i1=i1s,i1e
                icount=icount+1
                cvar1d(i1)=x_cntrl(icount)
                stddev1d(i1)=stddev(icount)
               end do 

          call perturb_1d (cv_name,i1s,i1e,stddev1d,correl,cvar1d,MPIRANK)

          deallocate(stddev1d)

               do i1=i1s,i1e
                icount_out=icount_out+1
                x(icount_out)=cvar1d(i1)
               end do 

             deallocate(cvar1d)
          else
          write(*,*) "IC DOES NOT HAVE APPROPRIATE ndim=",cvar_list(j)%ndim
          end if
         end if

!-- empirical model parameters
         if(cvar(i)%param) then
          cv_name=cvar(i)%name
          icount=icount+1
          cvarparm=x_cntrl(icount)
          icount_out=icount_out+1
         end if

!-- model bias
         if(cvar(i)%bias) then
         cv_name=cvar(i)%name
               do nbias=1,num_bias
           unitbfile=unitbias(nbias)
          if(cvar_list(j)%ndim.eq.4) then
          allocate(cvar4d(i1s:i1e,i2s:i2e,i3s:i3e,i4s:i4e))
!         allocate(stddev4d(i1s:i1e,i2s:i2e,i3s:i3e,i4s:i4e))

            do i4=i4s,i4e
             do i3=i3s,i3e
              do i2=i2s,i2e
               do i1=i1s,i1e
                icount=icount+1
                cvar4d(i1,i2,i3,i4)=x_cntrl(icount)
!               stddev4d(i1,i2,i3,i4)=stddev(icount)
               end do 
              end do 
             end do 
            end do 

!     call perturb_4d(cv_name,i1s,i1e,i2s,i2e,i3s,i3e,i4s,i4e, &
!                 stddev4d,correl,cvar4d,MPIRANK)
!         deallocate(stddev4d)
! For the moment we use default perturbations

            do i4=i4s,i4e
             do i3=i3s,i3e
              do i2=i2s,i2e
               do i1=i1s,i1e
                icount_out=icount_out+1
!               x(icount_out)=cvar4d(i1,i2,i3,i4)
               end do 
              end do 
             end do 
            end do 

          deallocate(cvar4d)
           else if(cvar_list(j)%ndim.eq.3) then
           allocate(cvar3d(i1s:i1e,i2s:i2e,i3s:i3e))
           allocate(stddev3d(i1s:i1e,i2s:i2e,i3s:i3e))

             do i3=i3s,i3e
              do i2=i2s,i2e
               do i1=i1s,i1e
                icount=icount+1
                cvar3d(i1,i2,i3)=x_cntrl(icount)
                stddev3d(i1,i2,i3)=stddev(icount)
               end do 
              end do 
             end do 

   call perturb_3d (cv_name,i1s,i1e,i2s,i2e,i3s,i3e,  &
                    stddev3d,correl,cvar3d,MPIRANK)

          deallocate(stddev3d)

             do i3=i3s,i3e
              do i2=i2s,i2e
               do i1=i1s,i1e
                icount_out=icount_out+1
                x(icount_out)=cvar3d(i1,i2,i3)
               end do 
              end do 
             end do 

           deallocate(cvar3d)
            else if(cvar_list(j)%ndim.eq.2) then
            allocate(cvar2d(i1s:i1e,i2s:i2e))
!           allocate(stddev2d(i1s:i1e,i2s:i2e))

              do i2=i2s,i2e
               do i1=i1s,i1e
                icount=icount+1
                cvar2d(i1,i2)=x_cntrl(icount)
!               stddev2d(i1,i2)=stddev(icount)
               end do 
              end do 

!  call perturb_2d (cv_name,i1s,i1e,i2s,i2e,stddev2d,correl,cvar2d,MPIRANK)

!         deallocate(stddev2d)

              do i2=i2s,i2e
               do i1=i1s,i1e
                icount_out=icount_out+1
!               x(icount_out)=cvar2d(i1,i2)
               end do 
              end do 

            deallocate(cvar2d)
             else if(cvar_list(j)%ndim.eq.1) then
             allocate(cvar1d(i1s:i1e))
             allocate(stddev1d(i1s:i1e))

               do i1=i1s,i1e
                icount=icount+1
                cvar1d(i1)=x_cntrl(icount)
                stddev1d(i1)=stddev(icount)
               end do 

   call perturb_1d (cv_name,i1s,i1e,stddev1d,correl,cvar1d,MPIRANK)

          deallocate(stddev1d)

               do i1=i1s,i1e
                icount_out=icount_out+1
                x(icount_out)=cvar1d(i1)
               end do 

             deallocate(cvar1d)
          else
          write(*,*) "BIAS DOES NOT HAVE APPROPRIATE ndim=",cvar_list(j)%ndim
          end if
               end do
       end if
!----------------------------------------------
        end if
        end do
       end do
!----------------------------------------------
!-- write current perturbation file (1d)
      WRITE(fileout,1000)  N
      CLOSE(iptrb)
      OPEN(iptrb,file=fileout,status='unknown',form='unformatted')
      WRITE(iptrb) idim_1d
      WRITE(iptrb) x
      CLOSE(iptrb,status='KEEP')

      END DO                              !!! loop over perturbations (ensembles)-N

 1000 format('x1d_',i4.4)
!--------------------------------------------

       if(MPIRANK.eq.0) then
       write(*,*) "end random"
       end if

!====================================================================
                                                                                               
#ifdef MPI_USE
      CALL MPI_FINALIZE(IERR)
#endif
                                                                                               
      stop
      end program random
!==========================================================================

                                subroutine perturb_1d                  &
   (cv_name,i1s,i1e,stddev1d,correl,cvar1d,MPIRANK)
                                                                                
!===  PROGRAM DOCUMENTATION BLOCK ============================= !                .      .    .                                       .
! PROGRAM: random_1d  APPLY KPZ equation to create Lyapunov-like vectors
! !          STEP 1:    INTERPOLATE SWM FIELDS TO LON-LAT GRADS FIELDS
!          STEP 2:    APPLY KPZ equation IN LON-LAT GRID
!          STEP 3:    INTERPOLATE LON-LAT TO SWM FIELDS
!
!   PRGMMR: ZUPANSKI        ORG: CIRA/CSU   DATE: 2004-06-25
!
! ABSTRACT: SMOOTH FORECAST PERTURBATIONS (FROM COLD START)
!
! PROGRAM HISTORY LOG:
! 2004-06-25  M.ZUPANSKI - 2-D shallow-water model scalars (no vectors yet)
! 2004-09-10  M.ZUPANSKI - Recursive Filter program
! 2004-09-18  M.ZUPANSKI - KPZ equation + random perturbations
!
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!   MACHINE:  IBM SP4
!===================================================================
      real,parameter    :: deg_incr=4.50   !!! average grid-point increment (degrees)
!----------------------------------------------------------------
!!    real,parameter    :: deg_incr=2.25   !!! average grid-point increment (degrees)
!----------------------------------------------------------------
      integer,parameter :: swm_in=22            !!! input SWM files  (1-D)
      real,parameter    :: PI=3.1415927
      real,parameter    :: DTR=3.1415927/180.
                                                                                
                                                                                
                                                                                
      real,parameter    :: qc_index=5.     !!! quality control index
      integer           :: i_lon,j_lat
      real              :: xlong
      real              :: latj
      integer           :: i_h
                                                                                
      integer :: i1s,i1e
      character*9 :: cv_name
      real,dimension(i1s:i1e) :: cvar1d
      real,dimension(i1s:i1e) :: stddev1d
      real :: correl
      real :: x_decorr,y_decorr
      integer :: MPIRANK

!--- interpolation weights and addresses
       real,dimension(:,:),allocatable      :: w1,w2,w3,w4
       integer,dimension(:,:),allocatable   :: k1,k2,k3,k4
       real,dimension(:),allocatable      :: z1,z2,z3,z4
       integer,dimension(:),allocatable   :: ic1,jc1
       integer,dimension(:),allocatable   :: ic2,jc2
       integer,dimension(:),allocatable   :: ic3,jc3
       integer,dimension(:),allocatable   :: ic4,jc4
                                                                                
!-- SWM variables
      real,dimension(:),allocatable ::   xlon,xlat
                                                                                
!-- LL variables
      real,dimension(:,:),allocatable :: x_LL_ptrb
      real,dimension(:),allocatable :: cvar1d_ptrb
                                                                                
!--- correlation variables
!---
      real,dimension(:,:),allocatable :: x_ptrb_x,x_ptrb_y
      real :: sumx,sumy
                                                                                
!======================================================
      x_decorr=correl
      y_decorr=correl
!----------------------------------------------------------
      i_swm=i1e-i1s+1      !!! No. of grid points
       if(MPIRANK.eq.0) then
      write(*,*) "START RANDOM_1d: cv_name=",cv_name," No OF PTS=",i_swm
      write(*,*) "deg_incr=",deg_incr
      write(*,*) "cvar1d min,max=",minval(cvar1d),maxval(cvar1d)
       end if
!----------------------------------------------------------
      i_lon=INT(360./deg_incr + 0.001)
      j_lat=INT(180./deg_incr + 0.001) + 1

       if(MPIRANK.eq.0) then
      write(*,*) "GRADS (lon-lat) grid dimensions: i,j=",i_lon,j_lat
      write(*,*) "GRADS (lon-lat) grid increment (degrees) =",deg_incr
       end if
                                                                                
!--- get input SWM variables (1-D)
                                                                                
      allocate(xlon(1:i_swm))          ! longitude (rad)
      allocate(xlat(1:i_swm))          ! latitude (rad)
                                                                                
         if(cv_name.eq.'height   ') then

      REWIND swm_in
        READ(swm_in)                !! h
        READ(swm_in)                !! u
        READ(swm_in)                !! v
        READ(swm_in)                !! psi
        READ(swm_in)                !! chi
        READ(swm_in)                !! rel
        READ(swm_in)                !! div
        READ(swm_in)                !! ke
        READ(swm_in)                !! eta
        READ(swm_in)                !! hs
        READ(swm_in) xlon           !! hlon
        READ(swm_in) xlat           !! hlat
        READ(swm_in)                !! ulon
        READ(swm_in)                !! ulat
        READ(swm_in)                !! urot
        READ(swm_in)                !! vrot

       if(MPIRANK.eq.0) then
       write(*,*) "INPUT: hlat min,max=",minval(xlat),maxval(xlat)
       write(*,*) "INPUT: hlon min,max=",minval(xlon),maxval(xlon)
       end if
                                                                                
       else

      REWIND swm_in
        READ(swm_in)                !! h
        READ(swm_in)                !! u
        READ(swm_in)                !! v
        READ(swm_in)                !! psi
        READ(swm_in)                !! chi
        READ(swm_in)                !! rel
        READ(swm_in)                !! div
        READ(swm_in)                !! ke
        READ(swm_in)                !! eta
        READ(swm_in)                !! hs
        READ(swm_in)                !! hlon
        READ(swm_in)                !! hlat
        READ(swm_in) xlon           !! ulon
        READ(swm_in) xlat           !! ulat
        READ(swm_in)                !! urot
        READ(swm_in)                !! vrot

       if(MPIRANK.eq.0) then
       write(*,*) "INPUT: ulat min,max=",minval(xlat),maxval(xlat)
       write(*,*) "INPUT: ulon min,max=",minval(xlon),maxval(xlon)
       end if
                                                                                
       end if

!-- Calculate (horizontal) interpolation weights
                                                                                
       allocate(w1(1:i_lon,1:j_lat))
       allocate(w2(1:i_lon,1:j_lat))
       allocate(w3(1:i_lon,1:j_lat))
       allocate(w4(1:i_lon,1:j_lat))
                                                                                
       allocate(k1(1:i_lon,1:j_lat))
       allocate(k2(1:i_lon,1:j_lat))
       allocate(k3(1:i_lon,1:j_lat))
       allocate(k4(1:i_lon,1:j_lat))
                                                                                
       w1=0.0
       w2=0.0
       w3=0.0
       w4=0.0

       k1=0
       k2=0
       k3=0
       k4=0
                                                                                
                                   call weights            &
       (i_lon,j_lat,i_swm,deg_incr,xlon,xlat,w1,w2,w3,w4,k1,k2,k3,k4)
                                                                                
       allocate(z1(1:i_swm))
       allocate(z2(1:i_swm))
       allocate(z3(1:i_swm))
       allocate(z4(1:i_swm))

       allocate(ic1(1:i_swm))
       allocate(ic2(1:i_swm))
       allocate(ic3(1:i_swm))
       allocate(ic4(1:i_swm))
       allocate(jc1(1:i_swm))
       allocate(jc2(1:i_swm))
       allocate(jc3(1:i_swm))
       allocate(jc4(1:i_swm))
                                                                                
       z1=0.0
       z2=0.0
       z3=0.0
       z4=0.0

       ic1=0
       ic2=0
       ic3=0
       ic4=0
       jc1=0
       jc2=0
       jc3=0
       jc4=0
                                                                                
                                   call weights_SWM        &
(i_lon,j_lat,i_swm,deg_incr,xlon,xlat,z1,z2,z3,z4,ic1,jc1,ic2,jc2,ic3,jc3,ic4,jc4)
                                                                                
       deallocate(xlon) ; deallocate(xlat)
                                                                                
!================================================================

!================================================================
!-- Interpolation and filtering (smoothing)
!---------------------------------------------------------------------------------
! (1) -- interpolate from SWM to LL grid ----------------
! call intpl_SWM_to_LL    (i_lon,j_lat,i_swm,w1,w2,w3,w4,k1,k2,k3,k4,cvar1d,x_LL)
                                                                                
      allocate(x_LL_ptrb(1:i_lon,1:j_lat))
      x_LL_ptrb(:,:)=0.0
                                                                                
! (2) -- Kardar-Parisi-Zhang (KPZ) equation ------------
                                                                                
                            call KPZ                          &
    (i_lon,j_lat,x_LL_ptrb,x_decorr,y_decorr,qc_index,MPIRANK)

               xfactor=qc_index/maxval(abs(x_LL_ptrb))
                 do j=1,j_lat
                  do i=1,i_lon
                   x_LL_ptrb(i,j)=xfactor*x_LL_ptrb(i,j)
                  end do
                 end do
                                                                                
                                                                                
! (3) -- Smoothing (filtering) of KPZ equation perturbations -----------------
                                                                                
  call recursive          (i_lon,j_lat,x_LL_ptrb,deg_incr,x_decorr,y_decorr,MPIRANK)

!-- "quality control" of created perturbations -----------------
               xfactor=qc_index/maxval(abs(x_LL_ptrb))
                 do j=1,j_lat
                  do i=1,i_lon
                   x_LL_ptrb(i,j)=xfactor*x_LL_ptrb(i,j)
                  end do
                 end do
                                                                                
!=========================================
      allocate(cvar1d_ptrb(i1s:i1e))
      cvar1d_ptrb(:)=0.0
                                                                                
! (4) -- interpolate from LL to SWM grid ----------------
                   call intpl_LL_to_SWM             &
(i_lon,j_lat,i_swm,z1,z2,z3,z4,ic1,jc1,ic2,jc2,ic3,jc3,ic4,jc4,x_LL_ptrb,cvar1d_ptrb)
!---------------------------------------------------------------------------------

 deallocate(w1) ; deallocate(w2) ; deallocate(w3) ; deallocate(w4)
 deallocate(k1) ; deallocate(k2) ; deallocate(k3) ; deallocate(k4)
 deallocate(z1) ; deallocate(z2) ; deallocate(z3) ; deallocate(z4)
 deallocate(ic1) ; deallocate(ic2) ; deallocate(ic3) ; deallocate(ic4)
 deallocate(jc1) ; deallocate(jc2) ; deallocate(jc3) ; deallocate(jc4)
                                                                                
     deallocate(x_LL_ptrb)
!=========================================
!-- perturbed control variable ------------

    do i=i1s,i1e
    cvar1d(i) = cvar1d(i) + stddev1d(i)*cvar1d_ptrb(i)
    end do
                                                                                
       if(MPIRANK.eq.0) then
      write(*,*) "cvar1d_ptrb min,max=",minval(cvar1d_ptrb),maxval(cvar1d_ptrb)
      write(*,*) "cvar1d min,max=",minval(cvar1d),maxval(cvar1d)
       end if
!---------------------------------------------------
     deallocate(cvar1d_ptrb)
                                                                                
     if(MPIRANK.eq.0) then
     write(*,*) " end perturb_1d for ",cv_name
     end if
!-----------------------------------------------------------------------
                            end subroutine perturb_1d
                                                                                
!============================================================

                                subroutine perturb_3d                  &
   (cv_name,i1s,i1e,i2s,i2e,i3s,i3e,stddev3d,correl,cvar3d,MPIRANK)

!===  PROGRAM DOCUMENTATION BLOCK =============================
!                .      .    .                                       .
! PROGRAM: perturb_3d APPLY KPZ equation to create Lyapunov-like vectors
!
!   PRGMMR: D. ZUPANSKI        ORG: CIRA/CSU   DATE: 2004-10-12
!
! ABSTRACT: SMOOTH FORECAST PERTURBATIONS (FROM COLD START)
!
! PROGRAM HISTORY LOG:
! 2004-06-25  M.ZUPANSKI - 2-D shallow-water model scalars (no vectors yet)
! 2004-09-10  M.ZUPANSKI - Recursive Filter program
! 2004-09-18  M.ZUPANSKI - KPZ equation + random perturbations
! 2004-10-12  D.ZUPANSKI - KPZ equation + random perturbations for 3d variables
!
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!   MACHINE:  IBM SP4
!===================================================================
      real,parameter    :: qc_index=5.     !!! quality control index

      integer :: i1s,i1e,i2s,i2e,i3s,i3e
      integer :: i,j,k,num_points
      character*9 :: cv_name
      real,dimension(i1s:i1e,i2s:i2e,i3s:i3e) :: cvar3d
      real,dimension(i1s:i1e,i2s:i2e,i3s:i3e) :: stddev3d
      real :: correl,x_decorr,y_decorr,z_decorr
      integer :: MPIRANK

      real,dimension(:,:,:),allocatable :: cvar3d_ptrb

      real :: DX,DY,DZ

!======================================================
      DX=60.
      DY=60.
      DZ=0.1
      x_decorr=correl
      y_decorr=correl
      z_decorr=0.
!----------------------------------------------------------
      num_points=(i1e-i1s+1)*(i2e-i2s+1)*(i3e-i3s+1)  !Tot numb of grid points
       if(MPIRANK.eq.0) then
      write(*,*) "START RANDOM_3d: cv_name=",cv_name," No OF PTS=",num_points
      write(*,*) "cvar3d min,max=",minval(cvar3d),maxval(cvar3d)
       end if
!----------------------------------------------------------
      allocate(cvar3d_ptrb(i1s:i1e,i2s:i2e,i3s:i3e))
      cvar3d_ptrb(:,:,:)=0.0

! (2) -- Kardar-Parisi-Zhang (KPZ) equation ------------
                            call KPZ_3d                       &
     (i1s,i1e,i2s,i2e,i3s,i3e,cvar3d_ptrb,x_decorr,y_decorr,z_decorr &
     ,qc_index,MPIRANK)

               xfactor=qc_index/maxval(abs(cvar3d_ptrb))
                do k=i3s,i3e
                 do j=i2s,i2e
                  do i=i1s,i1e
                   cvar3d_ptrb(i,j,k)=xfactor*cvar3d_ptrb(i,j,k)
                  end do
                 end do
                end do

       if(MPIRANK.eq.0) then
    write(*,*)"After KPZ: cvar3d_ptrb min,max=", &
    minval(cvar3d_ptrb),maxval(cvar3d_ptrb)
       end if
! (3) -- Smoothing (filtering) of KPZ equation perturbations -----------------

      call recursive_3d(i1s,i1e,i2s,i2e,i3s,i3e,cvar3d_ptrb, &
                        DX,DY,DZ,x_decorr,y_decorr,z_decorr,MPIRANK)

!-- "quality control" of created perturbations -----------------
               xfactor=qc_index/maxval(abs(cvar3d_ptrb))
                do k=i3s,i3e
                 do j=i2s,i2e
                  do i=i1s,i1e
                   cvar3d_ptrb(i,j,k)=xfactor*cvar3d_ptrb(i,j,k)
                  end do
                 end do
                end do

       if(MPIRANK.eq.0) then
       write(*,*) "After normalization: cvar3d_ptrb min,max=", &
       minval(cvar3d_ptrb),maxval(cvar3d_ptrb)
       end if
!=========================================

!-- perturbed control variable ------------


    do k=i3s,i3e
     do j=i2s,i2e
      do i=i1s,i1e
    cvar3d(i,j,k) = cvar3d(i,j,k) + stddev3d(i,j,k)*cvar3d_ptrb(i,j,k)
      end do
     end do
!     write(30) cvar3d(:,:,k)
    end do

       if(MPIRANK.eq.0) then
      write(*,*) "stddev3d min,max=",minval(stddev3d),maxval(stddev3d)
      write(*,*) "cvar3d_ptrb min,max=",minval(cvar3d_ptrb),maxval(cvar3d_ptrb)
      write(*,*) "cvar3d min,max=",minval(cvar3d),maxval(cvar3d)
       end if
!---------------------------------------------------
     deallocate(cvar3d_ptrb)

     if(MPIRANK.eq.0) then
     write(*,*) " end perturb_3d for ",cv_name
     end if
!-----------------------------------------------------------------------
                            end subroutine perturb_3d
!============================================================

                   subroutine intpl_SWM_to_LL                     &
        (i_lon,j_lat,i_swm,w1,w2,w3,w4,k1,k2,k3,k4,x,x_grads)

       integer :: i_lon,j_lat,i_swm
       integer :: ic,jc
       real,dimension(1:i_lon,1:j_lat)   :: w1,w2,w3,w4
       integer,dimension(1:i_lon,1:j_lat)   :: k1,k2,k3,k4
       real,dimension(1:i_swm) :: x                     !! input
       real,dimension(1:i_lon,1:j_lat) :: x_grads       !! output
!---------------------------------------------------------

         do ic=1,i_lon
          do jc=1,j_lat

      x_grads(ic,jc) = w1(ic,jc)*x(k1(ic,jc)) + w2(ic,jc)*x(k2(ic,jc)) &
                     + w3(ic,jc)*x(k3(ic,jc)) + w4(ic,jc)*x(k4(ic,jc))

          end do
         end do

          end subroutine intpl_SWM_to_LL
!==============================================================
                   subroutine intpl_SWM_to_LL_adj     &
        (i_lon,j_lat,i_swm,w1,w2,w3,w4,k1,k2,k3,k4,x_grads,x)

!---   INPUT: x_grads (LL-grid)
!---  OUTPUT: x       (SWM-grid)

       integer :: i_lon,j_lat,i_swm
       integer :: ic,jc
       real,dimension(1:i_lon,1:j_lat)   :: w1,w2,w3,w4
       integer,dimension(1:i_lon,1:j_lat)   :: k1,k2,k3,k4
       real,dimension(1:i_swm) :: x                     !! input
       real,dimension(1:i_lon,1:j_lat) :: x_grads       !! output
!---------------------------------------------------------

         do ic=1,i_lon
          do jc=1,j_lat

      x(k1(ic,jc))=x(k1(ic,jc))+w1(ic,jc)*x_grads(ic,jc)
      x(k2(ic,jc))=x(k2(ic,jc))+w2(ic,jc)*x_grads(ic,jc)
      x(k3(ic,jc))=x(k3(ic,jc))+w3(ic,jc)*x_grads(ic,jc)
      x(k4(ic,jc))=x(k4(ic,jc))+w4(ic,jc)*x_grads(ic,jc)

          end do
         end do

          end subroutine intpl_SWM_to_LL_adj

!============================================================

                   subroutine intpl_LL_to_SWM   &
(i_lon,j_lat,i_swm,w1,w2,w3,w4,ic1,jc1,ic2,jc2,ic3,jc3,ic4,jc4,x_LL,x)

       integer :: i_lon,j_lat,i_swm
       integer :: ic,jc
       real,dimension(1:i_swm)   :: w1,w2,w3,w4
       integer,dimension(1:i_swm)   :: ic1,jc1
       integer,dimension(1:i_swm)   :: ic2,jc2
       integer,dimension(1:i_swm)   :: ic3,jc3
       integer,dimension(1:i_swm)   :: ic4,jc4
       real,dimension(1:i_lon,1:j_lat) :: x_LL          !! input
       real,dimension(1:i_swm) :: x                     !! output
!---------------------------------------------------------

         do i=1,i_swm

      x(i) = w1(i)*x_LL(ic1(i),jc1(i)) + w2(i)*x_LL(ic2(i),jc2(i)) &
           + w3(i)*x_LL(ic3(i),jc3(i)) + w4(i)*x_LL(ic4(i),jc4(i))

         end do

          end subroutine intpl_LL_to_SWM
!==============================================================

!==============================================================

                         subroutine weights  &
       (i_lon,j_lat,i_swm,deg_incr,lon,lat,w1,w2,w3,w4,k1,k2,k3,k4)

       integer :: i_lon,j_lat,i_swm
       real,dimension(1:i_lon,1:j_lat)    :: w1,w2,w3,w4
       integer,dimension(1:i_lon,1:j_lat) :: k1,k2,k3,k4
       real,dimension(1:i_swm)    :: lon,lat
       real,dimension(:),allocatable    :: dl
       integer,dimension(:),allocatable    :: kw
       real    :: deg_incr
       real    :: lonc,latc
       real    :: xxlon,dxlon,dxlat
       real    :: s1,s2,s3,s4
       real    :: s_tot
       integer :: ic,jc
       real    :: loni
       real,parameter    :: PI=3.1415927
       real,parameter    :: DTR=PI/180.
       real,parameter    :: rearth=6371.
!--------------------------------------------------

     allocate(dl(1:i_swm))
     allocate(kw(1:i_swm))

            do ic=1,i_lon
           do jc=1,j_lat

      lonc=float(ic-1)*deg_incr*DTR          !!! East to West
      latc=float(jc-j_lat/2-1)*deg_incr*DTR  !!! 0 => Equator , >0 => NH , <0 => SH


         do i=1,i_swm
             if(lon(i).lt.0.) then
              loni=lon(i)+2.*PI
               else
                loni=lon(i)
                 end if
             xxlon=abs(lonc-loni)
             dxlon=min(xxlon,360.-xxlon)
             dxlat=latc-lat(i)
          dl(i)=rearth*sqrt( ( cos(latc) * dxlon )**2 + dxlat**2 )
          kw(i)=i
         end do

         call ORDER (dl,kw,i_swm,0)

!--- indexes
       k1(ic,jc)=kw(1)
       k2(ic,jc)=kw(2)
       k3(ic,jc)=kw(3)
       k4(ic,jc)=kw(4)

!--- weights
       if(dl(1).eq.0.) then
        w1(ic,jc)=1.
         w2(ic,jc)=0.
          w3(ic,jc)=0.
           w4(ic,jc)=0.
       elseif(dl(2).eq.0.) then
        w1(ic,jc)=0.
         w2(ic,jc)=1.
          w3(ic,jc)=0.
           w4(ic,jc)=0.
       elseif(dl(3).eq.0.) then
        w1(ic,jc)=0.
         w2(ic,jc)=0.
          w3(ic,jc)=1.
           w4(ic,jc)=0.
       elseif(dl(4).eq.0.) then
        w1(ic,jc)=0.
         w2(ic,jc)=0.
          w3(ic,jc)=0.
           w4(ic,jc)=1.
       else
        s1=1./dl(1)**4
         s2=1./dl(2)**4
          s3=1./dl(3)**4
           s4=1./dl(4)**4
!       s1=1./dl(1)**2
!        s2=1./dl(2)**2
!         s3=1./dl(3)**2
!          s4=1./dl(4)**2
       s_tot=s1+s2+s3+s4
         w1(ic,jc)=s1/s_tot
          w2(ic,jc)=s2/s_tot
           w3(ic,jc)=s3/s_tot
            w4(ic,jc)=s4/s_tot
       end if

           end do     !! j_lat
            end do     !! i_lon

       deallocate(dl)
       deallocate(kw)

          end subroutine weights

!==============================================================

                         subroutine weights_SWM  &
(i_lon,j_lat,i_swm,deg_incr,lon,lat,w1,w2,w3,w4,ic1,jc1,ic2,jc2,ic3,jc3,ic4,jc4)

       integer :: i_lon,j_lat,i_swm
       real,dimension(1:i_swm)    :: w1,w2,w3,w4
       integer,dimension(1:i_swm) :: ic1,jc1
       integer,dimension(1:i_swm) :: ic2,jc2
       integer,dimension(1:i_swm) :: ic3,jc3
       integer,dimension(1:i_swm) :: ic4,jc4
       real,dimension(1:i_swm)    :: lon,lat
       real,dimension(:),allocatable    :: dl
       integer,dimension(:),allocatable    :: kw
       real    :: deg_incr
       real    :: lonc,latc
       real    :: xxlon,dxlon,dxlat
       real    :: s1,s2,s3,s4
       real    :: s_tot
       integer :: ic,jc
       real    :: loni
       real,parameter    :: PI=3.1415927
       real,parameter    :: DTR=PI/180.
       real,parameter    :: rearth=6371.
!--------------------------------------------------

     allocate(dl(1:i_lon*j_lat))
     allocate(kw(1:i_lon*j_lat))

         do i=1,i_swm

             if(lon(i).lt.0.) then
              loni=lon(i)+2.*PI
               else
                loni=lon(i)
                 end if

          ij=0
            do ic=1,i_lon
           do jc=1,j_lat
          ij=ij+1

      lonc=float(ic-1)*deg_incr*DTR          !!! East to West
      latc=float(jc-j_lat/2-1)*deg_incr*DTR  !!! 0 => Equator , >0 => NH , <0 => SH

             xxlon=abs(lonc-loni)
             dxlon=min(xxlon,360.-xxlon)
             dxlat=latc-lat(i)
          dl(ij)=rearth*sqrt( ( cos(latc) * dxlon )**2 + dxlat**2 )
          kw(ij)=ij
           end do     !! j_lat
            end do     !! i_lon

         call ORDER (dl,kw,i_lon*j_lat,0)

!--- indexes
          ij=0
            do ic=1,i_lon
           do jc=1,j_lat
          ij=ij+1
          if(kw(1).eq.ij) then
           ic1(i)=ic
            jc1(i)=jc
             elseif(kw(2).eq.ij) then
              ic2(i)=ic
               jc2(i)=jc
                elseif(kw(3).eq.ij) then
                 ic3(i)=ic
                  jc3(i)=jc
                   elseif(kw(4).eq.ij) then
                    ic4(i)=ic
                     jc4(i)=jc
                   end if
           
           end do     !! j_lat
            end do     !! i_lon

!--- weights
       if(dl(1).eq.0.) then
        w1(i)=1.
         w2(i)=0.
          w3(i)=0.
           w4(i)=0.
       elseif(dl(2).eq.0.) then
        w1(i)=0.
         w2(i)=1.
          w3(i)=0.
           w4(i)=0.
       elseif(dl(3).eq.0.) then
        w1(i)=0.
         w2(i)=0.
          w3(i)=1.
           w4(i)=0.
       elseif(dl(4).eq.0.) then
        w1(i)=0.
         w2(i)=0.
          w3(i)=0.
           w4(i)=1.
       else
        s1=1./dl(1)**4
         s2=1./dl(2)**4
          s3=1./dl(3)**4
           s4=1./dl(4)**4
!       s1=1./dl(1)**2
!        s2=1./dl(2)**2
!         s3=1./dl(3)**2
!          s4=1./dl(4)**2
       s_tot=s1+s2+s3+s4
         w1(i)=s1/s_tot
          w2(i)=s2/s_tot
           w3(i)=s3/s_tot
            w4(i)=s4/s_tot
       end if

         end do           !! i_swm

       deallocate(dl)
       deallocate(kw)

          end subroutine weights_SWM
!================================================
                    subroutine ORDER &
       (V,IV,MAX,IREV)
!---------------------------------------
!  ORDER input vector V in ascending (IREV=0) or descending (IREV<>0) order
!  INPUT : V  - real vector (1:MAX)
!  INPUT : IV - integer (1:MAX), indexes 
! OUTPUT: V - ordered values
! OUTPUT: IV - ordered indexes
!--------------------------------------------------------------
      REAL V(*)
      INTEGER IV(*)
      integer :: I,IOFSET,LIM,ISW,NH
!---------------------------
      DO 10 I=1,MAX
        IV(I) = I
   10 CONTINUE
      IOFSET = MAX/2
   20 CONTINUE
      LIM = MAX - IOFSET
   30 CONTINUE
      ISW = 0
      DO 40 I=1,LIM
        IF(V(I).GT.V(I+IOFSET)) THEN
          VT = V(I)
          V(I) = V(I+IOFSET)
          V(I+IOFSET) = VT
          IVT = IV(I)
          IV(I) = IV(I+IOFSET)
          IV(I+IOFSET) = IVT
          ISW = I
        ENDIF
   40 CONTINUE
      LIM = ISW - IOFSET
      IF(ISW.NE.0) GO TO 30
      IOFSET = IOFSET/2
      IF(IOFSET.GT.0) GO TO 20

      IF(IREV.NE.0) THEN
!        REVERSE SORT ORDER...
         NH = MAX/2
         DO I=1,NH
            ITEMP = IV(I)
            IV(I) = IV(MAX+1-I)
            IV(MAX+1-I) = ITEMP
            TEMP = V(I)
            V(I) = V(MAX+1-I)
            V(MAX+1-I) = TEMP
         ENDDO
      ENDIF
      END subroutine ORDER

!==========================================================================
                                                                                
                            subroutine KPZ              &
    (imax,jmax,x_ptrb,r10,r20,qc_index,MPIRANK)
                                                                                
!-----
!----- KPZ parameters
!-----
   integer                         :: imax,jmax
   integer                         :: MPIRANK
   real,dimension(1:imax,1:jmax)   :: x_ptrb
   real :: r10,r20
                                                                                
!---
   integer :: i,j,k,it
   integer                         :: time,timex
   real,parameter    :: deg_incr=4.5
   real,parameter    :: PI=3.1415927
   real,parameter    :: DTR=3.1415927/180.
   real              :: qc_index
!===========================================================================
!-- initialize x_ptrb
     x_ptrb(:,:)=1.0
                                                                                
     time=5
!--- KPZ along longitude -------
     do i=1,imax
     timex=time
     qc_limit=qc_index*float(timex)/float(time)
     call KPZ_1d (jmax,x_ptrb(i,:),.false.,time,qc_limit,MPIRANK)
     end do
                                                                                
     time=INT(r10/5.+0.001)
!--- KPZ along latitude --------
!-- Poles -----
     x_ptrb(:,1)=0.0
     x_ptrb(:,jmax)=0.0
!-- Other latitudes
     do j=2,jmax-1
     jj=j-jmax/2-1
     timex=INT( float(time) * cos(deg_incr*float(jj)*DTR) )
     timex=max(1,timex)
     qc_limit=qc_index* float(timex)/float(time)
     call KPZ_1d (imax,x_ptrb(:,j),.true.,timex,qc_limit,MPIRANK)
     end do
!-------------------------------
                                                                                
      END subroutine KPZ

!==========================================================================

                            subroutine KPZ_1d           &
    (imax,x_ptrb,cyclic_bc,time,qc_limit,MPIRANK)

!-----------------------------------------------------------------
   integer                         :: imax
   integer                         :: time
   integer                         :: MPIRANK
   real,dimension(1:imax)          :: x_ptrb
   real                            :: qc_limit
   logical                         :: cyclic_bc

!-----
!----- KPZ parameters
!-----
   integer,parameter :: idiff=1              ! buffer points used in KPZ equation
   real,parameter    :: dt=0.1               ! time step length (nondimensional units)
   real,parameter    :: dx=1000.0            ! grid distance (nondimensional units)
   real,parameter    :: eps=dt/(dx**2)       ! KPZ epsilon parameter
   real,parameter    :: xlam=1.0
 
!--- random numbers
       integer :: i,j,k,it
       real,dimension(:),allocatable :: random_x
       real,dimension(:),allocatable :: v

       real  :: xfactor
!===========================================================================

     allocate(random_x(1:imax))
     allocate(v(1-idiff:imax+idiff))
     v(:)=0.0

!---------------------------------------
!---- main time loop -------------

      do it=1,time
                                                                                               
 do i=1,imax
  v(i)=x_ptrb(i)
   end do
                                                                                               
 if(cyclic_bc) then
      do i=1,idiff
       v(1-i)=x_ptrb(imax-i+1)
        v(imax+i)=x_ptrb(i)
         end do
 else
     do i=1,idiff
      v(1-i)=v(1+i)
       v(imax+i)=v(imax-i)
        end do
 end if

!-----------------------------
      mode2=MOD(imax,2)
       Nsave=imax/2 + mode2
        do l=1,Nsave
         call Box_Muller_polar(y1,y2)
          ii_tot=0
           do while (abs(y1).gt.3..or.abs(y2).gt.3.)
            ii_tot=ii_tot+1
             if(ii_tot.gt.10) exit
             call Box_Muller_polar(yy1,yy2)
            if(abs(yy1).le.3.) then
           y1=yy1
          end if
         if(abs(yy2).le.3.) then
        y2=yy2
       end if
      end do
     random_x(l)=y1
    random_x(Nsave+l-mode2)=y2
   end do
 
      do i=1,imax
     x_ptrb(i)=v(i) &
   +eps*(v(i+1)+v(i-1)-2.*v(i))+dt*0.5*xlam*v(i)+dt*sqrt(xlam)*v(i)*random_x(i)
  end do

       end do   !! time loop
!-------------------------------

   do i=1,imax
    x_ptrb(i)=sign(x_ptrb(i),random_x(i))
     end do

                if(maxval(abs(x_ptrb)).gt.0.) then
               xfactor=qc_limit/maxval(abs(x_ptrb))
              else
             xfactor=1.0
            end if

!!!         if(xfactor.gt.1) then
             do i=1,imax
              x_ptrb(i)=xfactor*x_ptrb(i)
               end do
!!!             end if

     deallocate(random_x)
 
!====================================================
      END subroutine KPZ_1d
!============================================================

                            subroutine KPZ_3d        &
     (i1s,i1e,i2s,i2e,i3s,i3e,x_ptrb,x_decorr,y_decorr,z_decorr &
     ,qc_index,MPIRANK)
                                                                                
!-----
!----- KPZ parameters
!-----
   integer                         :: i1s,i1e,i2s,i2e,i3s,i3e
   integer                         :: MPIRANK
   real,dimension(i1s:i1e,i2s:i2e,i3s:i3e)   :: x_ptrb
   real :: x_decorr,y_decorr,z_decorr
                                                                                
!---
   integer :: i,j,k,it
   integer                         :: time
   real              :: qc_index
!===========================================================================
!-- initialize x_ptrb
     x_ptrb(:,:,:)=1.0

            do k=i3s,i3e
                                                                                
     time=5
!--- KPZ along i1 -------
     imax=i2e-i2s+1
     do i=i1s,i1e
     call KPZ_1d (imax,x_ptrb(i,:,k),.false.,time,qc_index,MPIRANK)
     end do
                                                                                
     time=INT(x_decorr/5.+0.001)
!--- KPZ along i2 --------
     imax=i1e-i1s+1
     do j=i2s,i2e
     call KPZ_1d (imax,x_ptrb(:,j,k),.false.,time,qc_index,MPIRANK)
     end do
!-------------------------------
            end do
                                                                                
      END subroutine KPZ_3d
!============================================================

                      subroutine recursive                      &
  (imax,jmax,x,deg_incr,r10,r20,MPIRANK)
!------------------------------------------------------
      real,parameter    :: PI=3.1415927
      real,parameter    :: DTR=3.1415927/180.
      real,parameter    :: rearth=6371.

      integer :: imax,jmax
      integer :: MPIRANK
      real,dimension(1:imax,1:jmax) :: x
      real :: deg_incr
      real :: r10,r20

                                                                                                   
!--- correlation variables
      real,dimension(:),allocatable :: DX
      real                          :: DY
      real :: DD_lambda,DD_phi
      real,dimension(:),allocatable :: r_x,r_y
      integer,dimension(:),allocatable :: KX
      integer                          :: KY

      real,dimension(:,:),allocatable :: y

!--------------------------------------------------------------------
!   (elements obtained using SOAR-like compactly supported correlations)

        DD_lambda=deg_incr*DTR
        DD_phi=deg_incr*DTR

        DY=rearth*DD_phi

 allocate(DX(1:jmax))
        do j=1,jmax
         jj=j-jmax/2-1
!!         DX(j)=max(0.,rearth*cos(float(jj)*DD_phi)*DD_lambda)
           DX(j)=rearth*DD_lambda
             end do
!-------------------------------------------------------------

     allocate(y(1:imax,1:jmax))
     allocate(KX(1:jmax))

 if(DY.GT.0.) then
  KY = r20/DY + 1 
   else
    KY=0
     endif
  allocate(r_y(0:KY))
   r2=float(KY)*DY+0.1
  CALL CORR1D (r2,DY,KY,r_y)

        jloc=0

        do j=1,jmax

 if(DX(j).GT.0.) then
  KX(j) = r10/DX(j) + 1
   else
    KX(j)=0
     endif
         allocate(r_x(0:KX(j)))
          r1=float(KX(j))*DX(j) +0.1
          CALL CORR1D (r1,DX(j),KX(j),r_x)

      do i=1,imax

   y(i,j)=0.0
   do jj=max(1,j-KY),min(j+KY,jmax)
    ind_y=IABS(jj-j)
     do ii=max(1,i-KX(j)),min(i+KX(j),imax)
      ind_x=IABS(ii-i)
       y(i,j)=y(i,j) + r_x(ind_x)*r_y(ind_y)*x(ii,jj)
        end do
         end do

      end do
      deallocate(r_x)
        end do

      deallocate(r_y)

!-------------------------------------------
        x(:,:)=y(:,:)

     deallocate(y)
     deallocate(KX)
!------------------------------------------------------
                      end subroutine recursive

!============================================================

                      subroutine recursive_3d                   &
  (i1s,i1e,i2s,i2e,i3s,i3e,x,DX,DY,DZ,r1,r2,r3,MPIRANK)
!------------------------------------------------------

        integer :: i1s,i1e,i2s,i2e,i3s,i3e
        integer :: MPIRANK
        real,dimension(i1s:i1e,i2s:i2e,i3s:i3e) :: x
        real                   :: DX,DY,DZ
        real                   :: r1,r2,r3

        real,dimension(:,:,:),allocatable :: y
        real,dimension(:),allocatable :: r_x,r_y,r_z
        integer                          :: KX,KY,KZ
!-------------------------------------------------------------

     allocate(y(i1s:i1e,i2s:i2e,i3s:i3e))

   if(DX.GT.0.) then
     KX = r1/DX + 1
   else
     KX = 0
   endif
   if(DY.GT.0.) then
     KY = r2/DY + 1
   else
     KY = 0
   endif
   if(DZ.GT.0.) then
     KZ = r3/DZ
   else
     KZ = 0
   endif


       allocate(r_x(0:KX))
       CALL CORR1D (r1,DX,KX,r_x)

       allocate(r_y(0:KY))
       CALL CORR1D (r2,DY,KY,r_y)

       allocate(r_z(0:KZ))
       CALL CORR1D (r3,DZ,KZ,r_z)

        jloc=0
         do k=i3s,i3e
        do j=i2s,i2e
      do i=i1s,i1e

   y(i,j,k)=0.0

   do kk=max(i3s,k-KZ),min(k+KZ,i3e)
    ind_z=IABS(kk-k)
    do jj=max(i2s,j-KY),min(j+KY,i2e)
     ind_y=IABS(jj-j)
     do ii=max(i1s,i-KX),min(i+KX,i1e)
      ind_x=IABS(ii-i)
       y(i,j,k)=y(i,j,k) + r_x(ind_x)*r_y(ind_y)*r_z(ind_z)*x(ii,jj,kk)
     end do
    end do
   end do
         end do
        end do
       end do

!-------------------------------------------
        x(:,:,:)=y(:,:,:)

     deallocate(y)
     deallocate(r_x)
     deallocate(r_y)
     deallocate(r_z)
!------------------------------------------------------
                      end subroutine recursive_3d
!============================================================

                            SUBROUTINE CORR1D &
 (R1,DX,KX,rx)

  REAL					:: R1
  REAL					:: DX
  INTEGER				:: KX
  REAL, DIMENSION(0:KX)              :: rx
!------------------------------------------------------
!    define parameters


            XC=R1/2.
            XL=R1/5.


  if(KX.eq.0) then
   rx(0)=1.0
    else
     CALL SOAR (KX  ,KX,rx,DX,XC,XL)    
      end if

 end subroutine CORR1D

!=================================================================

               subroutine SOAR &
            (KBAND,KDIM,R,DD,CC,CL)

!----------------------------------------------------------------
!   parameters for compactly supported SOAR correlation function 
!   (Cohn and Gaspari,1999)
!----------------------------------------------------------------
!
!   c0=(1-|r|/L)*exp(-|r|/L)
!
!   f(r,L,c)=[c0-exp( (|r|-2c)/L )] / [1-exp( -2c/L )]          0<=|r|<=c
!
!   f(r,L,c)=[(2c-|r|)/L]*exp( -|r|/L )] / [1-exp( -2c/L )]     c<=|r|<=2c
!
!--------------------------------------------------------
!
!   Given R (decorr length), define c and L 
!
!------------------------------------------------------------------
!-----------------------------------
!     CC = R_DECORR/2.
!     CL = R_DECORR/5.
!     DD = DX, DY, or DZ
!   KBAND= Topelitz band
!   KDIM = DIMENSION of correlation function
!   (normally KDIM=KBAND)
!      R = correlation function
!-----------------------------------
 INTEGER			:: KDIM,KBAND
 REAL				:: DD,CC,CL
 REAL , DIMENSION (0:KDIM)	:: R

 INTRINSIC EXP,MAX

!========================================
    R(:)=0.0
    x1=2.*CC/CL 

         do IND=0,KBAND
    rr=IND*DD
    x2=rr/CL
      if(rr.LE.CC) then
 c0=(1.+x2)*exp(-x2)
 R(IND)=(c0-exp(x2-x1))/(1.-exp(-x1))
      else
 x3=MAX(0.,x1-x2)
 R(IND)=x3*exp(-x2)/(1.-exp(-x1))
      endif
         end do
!========================================

    end subroutine SOAR