module m_cloud
!
      use m_kinds
      use m_interp_NR, only : find_name
      implicit none
!
      private
      public :: clean_cloud
      public :: set_cloud
      public :: cloud
!
      integer, save :: ncloud
      integer, save :: k_iter
      integer, save :: prof_indexes(10)
      real(rkind3), allocatable, save :: cloud_table(:,:)
!
!  This module has software for determining either (1) if an IR observation 
!  is cloud affected or (2) if a MW observation is unduely affected 
!  by the surface or by precipitation.  
!
!  For IR, the effect is determined by first determining a probability that 
!  the observation is cloud affected.  That probability is defined by a 
!  piece-wise linear function of the cloud fraction. Separate probabilities
!  are determined for high, medium, and low clouds, based on the corresponding
!  cloud fractions (f). 
!
!  The probability function P(f) is defined by three parameters (a, b, c)
!  a: if f <= a,  then P(f)=0 (the scene is clear of this level cloud)
!  c: if f >= c,  then P(f)=1 (the scene is cloud contaminated)
!  b defines 2 ranges: 
!    if  a < f < b , P(f)= 0.5 * (f-a)/(b-a)
!    if  b < f < c , P(f)= 0.5 + 0.5*(f-b)/(c-b)
!  
!  A random number 0<x<1 is drawn.  If x<P(f), then the scene is contaminated
!  The values of a,b,c are provided in a resource file table that is read in,
!  with values dependent on observation type and on cloud level. The table
!  also has values of sigma: a cloud of some level is assigned a cloud top
!  pressure of ps*sigma.
!
!  For AMSU, the 4 values of f used to compute corresponding P(f) are:
!    (1) fraction of grid area covered by sea ice)
!    (2) fraction of grid area covered by land vs. ocean)
!    (3) 3 hour accumulated convective precip at the surface (meters)
!    (4) 3 hour accumulated stratiform precip at the surface (meters)  
!
!  Initial software by Ronald Errico  February 2008
!
!  
      contains
!
!
!  x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x 
!
      subroutine clean_cloud 
      deallocate (cloud_table)
      end subroutine clean_cloud
!
!  x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x 
!
      subroutine set_cloud (idatetime,nfields2d,nlevs,lrc,rcfile,d_type, &
                            f_names,box_size)
!
      integer*8 :: idatetime
      integer :: nfields2d
      integer :: nlevs
      integer :: lrc
      character(len= *) :: rcfile
      character(len= *) :: d_type
      character(len= 4) :: f_names(nfields2d)
      real(rkind1)      :: box_size
!
      integer :: n
      integer :: ibox_size
      integer :: irandom
      integer*8 :: i_random(2) ! only first value actually used
      character(len=8) :: cdum(3)
      character(len=4) :: cname(10)
      character(len=4) :: d_type_read
!
! Open resource file
      open(unit=lrc,file =trim(rcfile), form='formatted')
      print *,' '
      print *,'Set cloud table'
      print ('(3a,i3)'),' input file=',trim(rcfile),   &
              ' opened on unit=',lrc
!
! Read cloud table 
! First space down through records until 1st desired record found     
      do n=1,100   
        read (lrc,'(a4)') d_type_read
        if (d_type_read == d_type(1:4)) exit
      enddo        
!
      read (lrc,'(3(1x,a8,i5))') cdum(1),ncloud,cdum(2),irandom, &
                          cdum(3),ibox_size
      print ('(3(1x,a8,i5))'),cdum(1),ncloud,cdum(2),irandom, &
                          cdum(3),ibox_size
!
      box_size=real(ibox_size)
      allocate (cloud_table(4,ncloud))
      read  (lrc,'(1x,a8)') cdum(1)
      print ('(a8)'),cdum(1)
      do n=1,ncloud
        read  (lrc,'(1x,a8,2x,a4,4f6.2)') cdum(1),cname(n),cloud_table(1:4,n)
        print ('(1x,a8,2x,a4,4f8.4)'),cdum(1),cname(n),cloud_table(1:4,n)
      enddo
      close (lrc)
!
! Associate names of fields in table with names read from NR file
      call find_name (nfields2d,f_names(1:nfields2d),'pres', &
                      prof_indexes(1),.true.)
      do n=1,ncloud
        call find_name (nfields2d,f_names(1:nfields2d),cname(n), &
                        prof_indexes(n+1),.true.)
      enddo
!
! Set seed for random number generator to datetime
      i_random(1)=idatetime+irandom
      call random_seed (put=i_random(1:1))
      print ('(a,i12,a,i12)'),' Seed for random number generator = ' &
             ,i_random(1),'  idatetime=',idatetime 
!
! Set number of iterations required to detemine lcldtop from pcldttop 
      k_iter=int(log(real(nlevs+1))/log(2.)) + 1
!
      print ('(a,a5)'),' Cloud table and indexes filled for ',d_type
!
      end subroutine set_cloud 
!
!
!  x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x 
!
      subroutine cloud (nlevs,nfields2d,prof_2d,plevs,lcldtop,pcldtop)
!
      integer :: lcldtop
      integer :: nfields2d
      integer :: nlevs
      real(rkind1) :: plevs(nlevs)
      real(rkind1) :: prof_2d(nfields2d)
      real(rkind1) :: pcldtop
!
      integer :: i, j
      integer :: k_below
      real(rkind1) :: d
      real(rkind1) :: prob_cloud
      real(rkind1) :: sigma
      real(rkind1) :: x
!
      sigma=one_k1
      do i=1,ncloud
        j=prof_indexes(i+1)
        if (prof_2d(j) >= cloud_table(3,i) ) then
          prob_cloud=one_k1
        elseif (prof_2d(j) <= cloud_table(1,i) ) then
          prob_cloud=zero_k1
        else
          d=prof_2d(j)-cloud_table(2,i)  
          if (d >= zero_k1 ) then
            prob_cloud=0.5*(one_k3+d/(cloud_table(3,i)-cloud_table(2,i)))
          else 
            prob_cloud=0.5*(one_k3+d/(cloud_table(2,i)-cloud_table(1,i)))
          endif
        endif               
        call random_number(x)  
        if (x < prob_cloud) then
          sigma=cloud_table(4,i)
        endif
        if (sigma /= one_k1) exit
      enddo
!     
      pcldtop=sigma*prof_2d(prof_indexes(1))
!
! Find the grid point index for the pressure level just above pcldtop by
! successively dividing the range it is in by factors of 2.
      lcldtop=1
      k_below=nlevs
      if (pcldtop >= plevs(nlevs)) then
        lcldtop=nlevs
      else
        do i=1,k_iter
          j=(lcldtop+k_below)/2 
          if (pcldtop >= plevs(j)) then
            lcldtop=j
          else
            k_below=j
          endif
        enddo
      endif     ! test if cloud top between surface and lowest level
!
      end subroutine cloud
! 
      end module m_cloud