MODULE aodvismod   
!====================================================================
!
! Name: aodvismod
!
! Function:
!	Evaluate Visibility
!
! Description:
!       Visibility is evaluated based on visible aerosol optical depth 
!       and planetary boundary layer depth
!       visibility due to fog will be added shortly into algorithm
!
! Reference:
!
! Inputs:
!	Input passed though GEOCAT structure: satellite, nwp,out2
!
! Outputs:
!
!	out2(Algo_Num)%vis_fromaod = surface visibility
!
! Public Routines: 
!	aodvis
!
! Public Variables:
!
! Dependencies: 
!	Geocat objects - visible aersol optical depth, 
!       nwp planetary boundary
!	layer depth, low cloud
!       to be added
!       ,fog probability, fog depth
!
! Restrictions:
!       No paths greater than 256 chars
!
! History:
! 7/14/2010 - Allen Lenzen - Version 1 - fog part of visibility to be added next 
!             week
!             when aod_vis is available in out2, algorithm will use that 
!             instead of using readaod routine
!
!====================================================================

!--- module use statements
use ALGORITHM_MODULE_USAGE

implicit none

private

!--- module-level subroutines and functions
public:: aodvis   ! main entry subroutine
                             
!--- module-level variables
real, allocatable :: glon(:,:),glat(:,:),gaod(:,:)

contains

SUBROUTINE aodvis(Algo_num)   
!====================================================================
!
! Name: aodvis
!
! Function:
!	Main visibility algoritm
!
! Description:
!	
!       Visibility is evaluated based on visible aerosol optical depth
!       and planetary boundary layer depth
!
! Reference:
!
! Calling Sequence: CALL aodvis(ialgo)
!
! Inputs: 
!	Input passed though GEOCAT structure: satellite, nwp,out2
!
! Outputs:
!	out2(Algo_Num)%vis_fromaod = surface visibility
!
! Dependencies: 
!	Input passed though GEOCAT structure: satellite, nwp,out2
!       bias corrections read in readvisbiascor
!
! Restrictions:
!
! History:
!
!====================================================================

  !--- arguments
  INTEGER(kind=int4), intent(in) :: Algo_num  ! algorithm index

  !--- local variables
  INTEGER(kind=int4) :: iline              ! line index
  INTEGER(kind=int4) :: ielem              ! element index
  INTEGER(kind=int4) :: NWP_Lon_Cell
  INTEGER(kind=int4) :: NWP_Lat_Cell
  REAL(kind=real4), dimension(:,:), POINTER :: Visibility,Aodvisible
  REAL(kind=real4), dimension(:,:), POINTER :: zpblpnt
  REAL(kind=real4) visib
  integer(kind=int4) :: visib_cat
  integer(kind=int1), dimension(:,:), POINTER :: pblcldpnt
  REAL(kind=real4) :: Pblh
  REAL (kind=real4),allocatable :: aaod_vis(:,:)
  logical notclear                         
! cloudy, probably cloudy, probably clear
  logical clear
! cloud mask 0-cloudy, 1-probably cloudy, 2-probably clear, 3-clear
! specifly number of months and catagories
  INTEGER(kind=int4),parameter :: nmonbias=12,ncatbias=4
  REAL (kind=real4),dimension(nmonbias,ncatbias),save :: modis_bias,modis_slope
  LOGICAL first
! paths to input data 
  character *200 dirbias,diraod
  integer il,im,imon
! so only read input once
  save first
  data first/.true./
! set up bias correction table on first call
  if(first)then
    print *,'do first',first
    first=.false.
    dirbias=TRIM(out2(Algo_Num)%ancil_subdir)
    diraod=TRIM(out2(Algo_Num)%ancil_subdir)
    write(6,*)'dirbias',dirbias
    call flush(6)
    call readvisbiascor (dirbias,nmonbias,ncatbias,modis_bias,modis_slope)
    do im=1,12
      do il=1,4
        print *,'modis_bias',im,il,modis_bias(im,il),modis_slope(im,il)
      end do
    end do
  endif
  imon=timestr%month
  !------------------------------
  ! Alias surface visibility
  Visibility => out2(Algo_Num)%vis_fromaod
  ! alias zpbl, pblcld, zpbl
  zpblpnt => out2(Algo_Num)%zpbl_vis
  pblcldpnt => out2(Algo_Num)%pblcld
  ! intialize to missing
  zpblpnt=missing_value_real4
  Visibility=missing_value_real4
  pblcldpnt=0
 
  !
  ! Alias aerosol optical depth - visible
  !   Aodvisible => out2(Algo_num)%aod_vis
  !   Aodvisible => out2(3)%aod_vis
  !   output input aod_vis for diagnostics for Version 1
  Aodvisible => out2(Algo_Num)%modisaod_abi
  Aodvisible=missing_value_real4
  ! diagnostic print for Version 1
  write(6,*)'Algo_num',Algo_num
  call flush(6)

  !--- preprocessing
  
  !--- main pixel loop
  allocate (aaod_vis(sat%nx,sat%ny))
  ! this is for Version 1 will use aod_vis when available in out2 in future
  call readaod(diraod,sat%xstart,sat%ystart,sat%nx,sat%ny,aaod_vis,sat%lat, &
    sat%lon)
  write(6,*)'nx',sat%nx,sat%ny
  write(6,*)'nx0',sat%nx0,sat%ny_all
  write(6,*)'xstart',sat%xstart,sat%ystart
  write(6,*)'iseg',sat%iseg
  call flush(6)
  line_loop: do iline=1, sat%ny
     element_loop: do ielem=1, sat%nx
        
        !--- check for space pixel
        if (sat%space_mask(ielem,iline) == sym%SPACE) then
           cycle
        endif

        !--- define aliases for this loop iteration
        
        !--- perform operations on the current pixel
        NWP_Lon_Cell = sat%x_nwp(ielem,iline)
        NWP_Lat_Cell = sat%y_nwp(ielem,iline)
        Pblh = nwp%dat(NWP_Lon_Cell,NWP_Lat_Cell)%zpbl

        !--- write output to the out2 struct
        Aodvisible(ielem,iline)=aaod_vis(ielem,iline)
        ! logical for notclear and clear
        notclear=sat%cldmask(ielem,iline).ne.3
        clear=sat%cldmask(ielem,iline).eq.3
        if (Aodvisible(ielem,iline) > 0.0 .and. clear) then
          Visibility(ielem,iline)=3.9/(Aodvisible(ielem,iline)/Pblh)
          !bias correction of visibility
          visib=Visibility(ielem,iline)/1000.
          visib_cat=0
          if(visib .ge. 30.)then
            visib_cat=1
          elseif(visib .ge. 10.)then
            visib_cat=2
          elseif(visib .ge. 2.)then
            visib_cat=3
          elseif(visib .gt. 0.)then
            visib_cat=4
          endif
          if(visib_cat .ne. 0 )then
            visib=.4*visib+.6*(modis_bias(imon,visib_cat) +  &
              modis_slope(imon,visib_cat)*visib)
          endif
        
          !   max visibility
          Visibility(ielem,iline)=min(30.e3,Visibility(ielem,iline))

        endif
        zpblpnt(ielem,iline)=nwp%dat(NWP_Lon_Cell,NWP_Lat_Cell)%zpbl
        ! set up flag for pblcld for Version 1 diagnostics
        if(sat%cldz(ielem,iline) > 0.0 )then
          if(sat%cldz(ielem,iline)-nwp%dat(NWP_Lon_Cell,NWP_lat_Cell)%zsfc &
            < zpblpnt(ielem,iline) ) then
          pblcldpnt(ielem,iline)=1
          endif
        endif

     end do element_loop
  end do line_loop
  deallocate (aaod_vis)

  !--- postprocessing

end SUBROUTINE aodvis
! end main subroutine
!
!====================================================================
!
! Name: Readvisbiascor
!
! Function:
!	Reads bias corrections for visibility
!
! Description:
!	Reads the coefficients into a properly-shaped 2D array
!
!
! Reference:
!
! Inputs:
!       directory to find bias correction coefficients
!       dimension of bias
!
! Outputs:
!
!	Bias correction coefficients (Bias and slope)
!
! Public Routines: 
!
! Public Variables:
!
! Dependencies: 
!	 
!	
!
! Restrictions:
!       No paths greater than 256 chars
!
! History:
SUBROUTINE readvisbiascor(dirbias,nmonbias,ncatbias,modis_bias,modis_slope)
implicit none
integer (kind=int4)nmonbias,ncatbias,nmonbiasi,ncatbiasi
real (kind=real4),dimension(nmonbias,ncatbias) :: modis_bias,modis_slope
character *200 filebias,dirbias*(*)
!filebias='/home/lenzen/geocat/vis/MODIS_ASOS_bias_correction.dat'
filebias=trim(dirbias)//'/MODIS_ASOS_bias_correction.dat'
open(10,file=filebias,form='unformatted')
! read in dimension of bias corrrections
read(10)nmonbiasi,ncatbiasi
print *,'nmonbias',nmonbiasi,ncatbiasi
! check dimensions
if(nmonbiasi.ne.nmonbias.or.ncatbiasi.ne.ncatbiasi)then
  print *,'error dimensions wrong readbiascor'
  stop
endif
! read bias and slope
  
read(10)modis_bias,modis_slope
close(10)
return
end SUBROUTINE readvisbiascor
! 
! --------------------------------
! this subroutine is for Version 1
! will use aod_vis when available in out2
SUBROUTINE readaod(diraod,xstart,ystart,nx,ny,aod,satlat,satlon)
implicit none
character *160 fileaod,yearjulday*7,diraod*(*)
integer nc,nr,nx,ny,ii,jj,i,j
real aod(nx,ny),satlat(nx,ny),satlon(nx,ny)
integer xstart,ystart
! so do read only once
logical doread
data doread/.true./
save doread,nc,nr
write(6,*)'xstart',xstart,ystart,'nx',nx,ny
call flush(6)
if(doread)then
  doread=.false.
  write(yearjulday,'(i4.4,i3.3)')timestr%year,timestr%jday
  fileaod=trim(diraod)//'/ABI_aot_'//yearjulday//'_conus.dat'

  open(10,file=fileaod,form='unformatted')
  read(10)nc,nr
  write(6,*)'readaod ncnr',nc,nr,'nx',nx,ny
  call flush(6)

  write(6,*)'readaod year ',timestr%year,' day ',timestr%jday
  write(6,*)'readaod yearjulday',yearjulday
  call flush(6)
  allocate (glon(nc,nr),glat(nc,nr),gaod(nc,nr))
  ! read in lon, lat and aod
  read(10)glon,glat,gaod
  close(10)
endif
! mv current aod to aod array
do j=1,ny
  jj=ystart+j-1
  do i=1,nx
    ii=xstart+i-1
    if(satlat(i,j).ne.-99999.)then
    if(abs(satlat(i,j)-glat(ii,jj))>.1)then
      write(6,*)'lat not equal',i,j,satlat(i,j),glat(ii,jj),'ii',ii,jj, &
        'ystart',ystart
    endif
    if(abs(satlon(i,j)-glon(ii,jj))>.1)then
       write(6,*)'lon not equal',i,j,satlon(i,j),glon(ii,jj),'ii',ii,jj, &
         'ystart',ystart
    endif
    endif
    aod(i,j)=gaod(ii,jj)
  end do
end do
return
end SUBROUTINE readaod
! end readaod subroutine


end MODULE aodvismod 

