MODULE visibmod   
!====================================================================
!
! Name: visibmod
!
! Function:
!	Evaluate Visibility
!
! Description:
!       In clear areas, aerosol Visibility is evaluated based on visible 
!       aerosol optical depth and planetary boundary layer depth.
!       In cloudy areas, fog visibility is evaluated based on cloud optical 
!       depth and fog probability and fog depth.
!       Cloud optical depth and effective cloud radius will be read from 
!       MODIS_MOD06 files.
!       Aerosol optical depth will be read from MODIS_MOD04 files and then
!       interpolated to pixel resolution.
!
! Reference:
!
! Inputs:
!	Input passed though GEOCAT structure: satellite, nwp, out2
!       Cloud optical depth and effective cloud radius will be read from
!        MODIS_MOD06 files.
!       Aerosol optical depth will be read from MODIS_MOD04 files.
!
! Outputs:
!
! 	out2(Algo_Num)%visib = surface visibility
!
! Diagnostic outputs:
!	out2(Algo_Num)%vis_type = surface visibility type
!                             0 - visibility is from aerosol optical depth
!                             1 - visibility is from fog visibility estimate
!
!	out2(Algo_Num)%Fguess_fog_visib = first guess visibility from fog
!	out2(Algo_Num)%fog_visib = bias corrected visibility from fog
!	out2(Algo_Num)%Fguess_aod_visib = first guess visibility from aod
!	out2(Algo_Num)%aod_visib = bias corrected visibility from aod 
!	out2(Algo_Num)%pblcld = 1 from cloud with planetary boundary layer
!       out2(Algo_Num)%zpbl_vis - Planetary boundary layer height
!
! Public Routines: 
!	visib_main
!
! Public Variables:
!
! Dependencies: 
!	Geocat objects - visible aersol optical depth, 
!       nwp planetary boundary layer depth,
!       fog probability, fog depth
!       cloud mask
!       for pblcld flag need cloud top height,surface height
!       cloud optical thickness
!       effective cloud drop radius
!
! Restrictions:
!       No paths greater than 256 chars
!       at present image not over a pole or over date line
!
! History:
! 7/14/2010 - Allen Lenzen - Version 1 - fog part of visibility to be added 
! august 2010 - added fog visibility
! september 2010 -added visible aerosol optical depth from GOES East aerosol
! /smoke product (GASP) ( now removed)
!  added readgaspaod subroutine (now removed)
!             when aod_vis is available in out2, algorithm will use that 
!             instead of using readgaspaod routine (Now removed)
! April 2011 - add code to read cloud optical thickness (COT) from MODIS_MOD06
!              files
!            - add code to read cloud effective radius (reff) from MODIS_MOD06
!              files
!            - add code to read aerosol optical depth (AOD) from MODIS_MOD04
!              files 
! May 2011   - add interpolation for aerosol optical depth from MOD04 resolution
!              to image resolution
!            - added updated look up tables (LUT) for aerosol visibility based
!              on AOD value and
!              fog visibility based on COT value
!            - change fog probability to 80 percent for fog present 
!            - aod interpolation changed to handle Aqua data properly
!==================================================================== 
!--- module use statements
use ALGORITHM_MODULE_USAGE
use ATMOSPHERE

implicit none

private

!--- module-level subroutines and functions
public:: visib_main   ! main entry subroutine
                             
!--- module-level variables

contains

SUBROUTINE visib_main(Algo_num)   
!====================================================================
!
! Name: visib
!
! Function:
!	Main visibility algoritm
!
! Description:
!	
!       In clear areas, aerosol Visibility is evaluated based on visible
!        aerosol optical depth and planetary boundary layer depth.
!       In cloudy areas, fog visibility is evaluated based on cloud optical 
!       depth and fog probability and fog depth.
!
! Reference:
!
! Calling Sequence: CALL vis(ialgo)
!
! Inputs: 
!	Input passed though GEOCAT structure: satellite, nwp, out2
!
! Outputs:
!	out2(Algo_Num)%visib = surface visibility

! Diagnostic outputs:
!	out2(Algo_Num)%vis_type = surface visibility type
!                             0 - visibility is from aerosol optical depth
!                             1 - visibility is from fog visibility estimate
!
!	out2(Algo_Num)%Fguess_fog_visib = first guess visibility from fog
!	out2(Algo_Num)%fog_visib = bias correcte visibility from fog
!	out2(Algo_Num)%Fguess_aod_visib = first guess visibility from aod
!	out2(Algo_Num)%aod_visib = bias corrected visibility from aod 
!	out2(Algo_Num)%pblcld = 1 from cloud with planetary boundary layer
!       out2(Algo_Num)%zpbl_vis - Planetary boundary layer height
!                                             
!
! Dependencies: 
!	Input passed though GEOCAT structure: satellite, nwp,out2
!	Geocat objects - visible aersol optical depth, 
!       nwp planetary boundary layer depth,
!       fog probability, fog depth
!       cloud mask
!       for pblcld flag need cloud top height,surface height
!       bias corrections read in with subroutine readvisbiascor
!
! Restrictions:
!
! History:
!
!====================================================================
use ATMOSPHERE

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

  !--- local variables
  INTEGER(kind=int4) :: Stat               ! status
  INTEGER(kind=int4) :: iline              ! line index
  INTEGER(kind=int4) :: ielem              ! element index
  INTEGER(kind=int4) :: NWP_Lon_Cell
  INTEGER(kind=int4) :: NWP_Lat_Cell
  integer(kind=int4) idim(2)
  REAL(kind=real4), dimension(:,:), POINTER :: Visibility
  REAL(kind=real4), dimension(:,:), pointer :: aod_visible
  REAL(kind=real4), dimension(:,:), POINTER :: Fguess_aod_vis
  REAL(kind=real4), dimension(:,:), POINTER :: Fguess_fog_vis
  REAL(kind=real4), dimension(:,:), POINTER :: fog_visib
  REAL(kind=real4), dimension(:,:), POINTER :: aerosol_visib
  REAL(kind=real4), dimension(:,:), POINTER :: zpblpnt
  REAL(kind=real8), dimension(:,:), POINTER,save :: COT
  REAL(kind=real8), dimension(:,:), POINTER,save :: reff
  REAL(kind=real8), dimension(:,:), POINTER,save :: CWP 
  REAL(kind=real8), dimension(:,:), POINTER,save :: AOD04
  REAL(kind=real8), dimension(:,:), POINTER,save :: MOD06_Lon 
  REAL(kind=real8), dimension(:,:), POINTER,save :: MOD06_Lat
  REAL(kind=real8), dimension(:,:), POINTER,save :: MOD04_Lon
  REAL(kind=real8), dimension(:,:), POINTER,save :: MOD04_Lat
  character(Len=256) :: Err_Message
  integer (kind=int1), dimension(:,:), POINTER :: Visib_type
  REAL(kind=real4) visib,visib_fog,aodvis
  integer(kind=int4) :: aod_cat,visib_cat
  integer(kind=int1), dimension(:,:), POINTER :: pblcldpnt
  REAL(kind=real4) :: Pblh
  logical notclear                         
! cloudy, probably cloudy, probably clear
  logical clear
  logical modis
  logical abi
! 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
! aod and fog bias corrections and slope
  REAL (kind=real4),dimension(nmonbias,ncatbias),save :: aod_bias
  REAL (kind=real4),dimension(nmonbias,ncatbias),save :: aod_slope
  REAL (kind=real4),dimension(nmonbias,ncatbias),save :: fog_bias
  REAL (kind=real4),dimension(nmonbias,ncatbias),save :: fog_slope
  integer(kind=int4) modissource6,modissource4
  LOGICAL first
! paths to input data 
  character *200 dirbias
  integer(kind=int4) imon
  integer (kind=int4)nc,nr,nc6,nr6
  integer (kind=int4) nc4,nr4
! so only read input once
  save first,modis,abi
  data first/.true./
  real(kind=real8), allocatable :: lwp(:,:),fog_depth(:,:)
  real(kind=real8) :: h2o_density
  real(kind=real4) :: lwc
  real(kind=real4) :: maxvisib
  data maxvisib/60./
!  data maxvisib/30./
  integer(kind=int4) ifog ! algorithm for fog
  integer(kind=int4) :: Lines_Add
! this needs to be changed if fog is not the fourth algorithm
  ifog=4
  allocate (lwp(sat%nx,sat%ny),fog_depth(sat%nx,sat%ny))
  lwp=missing_value_real8
! set up bias correction table on first call
  if(first)then
!    print *,'scinfo',scinfo(sc_ind)%name
    first=.false.
    dirbias=TRIM(out2(Algo_Num)%ancil_subdir)
    if(scinfo(sc_ind)%series.eq.'modis')then
     modis=.true.
     abi=.false.
    elseif(scinfo(sc_ind)%series.eq.'goesr')then
     modis=.false.
     abi=.true.
    else
      print *,'not valid series',scinfo(sc_ind)%series
      stop
    endif
    write(6,*)'VISIBILITY read bias, slope files (INFORMATION) maxvisib',maxvisib
    call readvisbiascor (dirbias,nmonbias,ncatbias,aod_bias,aod_slope, &
       fog_bias,fog_slope,modis,abi)
  endif
  imon=timestr%month
  !------------------------------
  ! Alias surface visibility
  Visibility => out2(Algo_Num)%visib
  Visib_type => out2(Algo_Num)%vis_type
  ! 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
  Visib_type=missing_value_int1
  pblcldpnt=0

! alias First guess visibilities

  Fguess_aod_vis => out2(Algo_Num)%fguess_aerosol_visib
  Fguess_fog_vis => out2(Algo_Num)%fguess_fog_visib
! alias visibilities
  fog_visib => out2(Algo_Num)%fog_visib
  aerosol_visib => out2(Algo_Num)%aerosol_visib
  aerosol_visib=missing_value_real4
  Fguess_fog_vis=missing_value_real4
  Fguess_aod_vis=missing_value_real4
  fog_visib=missing_value_real4
 
  !
  ! Alias aerosol optical depth - visible
  Aod_visible => out2(Algo_num)%aod_visible
  aod_visible=missing_value_real4

  !--- preprocessing
  
  !--- main pixel loop
  ! this is for Version 1 will use aod_vis when available in out2 in future
  nc=sat%nx
  nr=sat%ny
! in future for ABI
  if(abi)then
!   AOD, ((COT,reff) or LWP ) from sat or out2 structure
!   fog_probability, fog_depth from sat or out2 structure
  endif
  
! MODIS NOW
if(modis)then
! read in COT from MODIS_MOD06

  if(scinfo(sc_ind)%name.eq.'Terra')then
    modissource6=SOURCE_MODIS_MOD06
    modissource4=SOURCE_MODIS_MOD04
  else
    modissource6=SOURCE_MODIS_MYD06
    modissource4=SOURCE_MODIS_MYD04
  endif
  CALL Get_Atmosphere_Data(modissource6, &
       'Cloud_Optical_Thickness', &
       COT, &
       Native_Lon=MOD06_Lon, &
       Native_Lat=MOD06_Lat, &
       Stat=Stat)

  IF( Stat /= STATUS_OK ) THEN 
     WRITE (Err_Message, *) &
          "COT: Error from Get_Atmosphere_Data: ", &
          Stat
     CALL Display_Message('test_modis_ingest', Err_Message, Sym%FAILURE)
     RETURN
  END IF
  

! read in reff from MODIS_MOD06
  CALL Get_Atmosphere_Data(modissource6, &
       'Cloud_Effective_Radius', &
       reff, &
       Native_Lon=MOD06_Lon, &
       Native_Lat=MOD06_Lat, &
       Stat=Stat)

  IF( Stat /= STATUS_OK ) THEN 
     WRITE (Err_Message, *) &
          "COT: Error from Get_Atmosphere_Data: ", &
          Stat
     CALL Display_Message('test_modis_ingest', Err_Message, Sym%FAILURE)
     RETURN
  END IF
! read in cloud water path  from MODIS_MOD06
  CALL Get_Atmosphere_Data(modissource6, &
       'Cloud_Water_Path', &
       CWP, &
       Native_Lon=MOD06_Lon, &
       Native_Lat=MOD06_Lat, &
       Stat=Stat)
  idim=shape(MOD06_Lon)
  nc6=idim(1)
  nr6=idim(2)
  

  IF( Stat /= STATUS_OK ) THEN 
     WRITE (Err_Message, *) &
          "COT: Error from Get_Atmosphere_Data: ", &
          Stat
     CALL Display_Message('test_modis_ingest', Err_Message, Sym%FAILURE)
     RETURN
  END IF
! read in AOD from MODIS_MOD04
  !--------------------------------------------------
  ! AOD at native res
  !--------------------------------------------------
  
  Lines_Add=2
! this could probably be 1 but 2 is safer and very slight effect 
! on speed or memory size
  CALL Get_Atmosphere_Data(modissource4, &
       'Optical_Depth_Land_And_Ocean', &
       AOD04, &
        Pad_Lines=Lines_Add, &
       Native_Lon=MOD04_Lon, &
       Native_Lat=MOD04_Lat, &
       Stat=Stat)
  idim=shape(MOD04_lon)
  nc4=idim(1)
  nr4=idim(2)
  call aod04to06(nc6,nr6,MOD06_Lon,MOD06_Lat,&
  nc4,nr4,MOD04_Lon,MOD04_Lat,nc,nr,aod_visible,AOD04,sat%ystart,sat%ny_all)

  IF( Stat /= STATUS_OK ) THEN 
     WRITE (Err_Message, *) &
          "AOD small: Error from Get_Atmosphere_Data: ", &
          Stat
     CALL Display_Message('test_modis_ingest', Err_Message, Sym%FAILURE)
     RETURN
  END IF
endif ! MODIS

 
  h2o_density=998.2071 ! kg/m3 at 1 atm/20 degrees C 
  lwc=0.3 ! coefficient for fogdepth from lwp
  fog_depth=missing_value_real4
  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
        ! calculate lwp (liquid water path)
        if(COT(ielem,iline)>0 .and. reff(ielem,iline)>0)then
          lwp(ielem,iline)=(5./9.)*h2o_density*COT(ielem,iline)* &
            reff(ielem,iline)*1.e-6*1.e3
          fog_depth(ielem,iline)=lwp(ielem,iline)/lwc ! units meters
        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
        ! logical for notclear and clear
        notclear=sat%cldmask(ielem,iline).ne.3
        clear=sat%cldmask(ielem,iline).eq.3
        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
        if (Aod_visible(ielem,iline) > 0.0 .and. clear) then
!         use WMO value now which is 3.0
          Visibility(ielem,iline)=3./(Aod_visible(ielem,iline)/(Pblh/1000.))
          !bias correction of visibility
          ! now visib in killometers
          visib=Visibility(ielem,iline)
          Fguess_aod_vis(ielem,iline)=visib
          aodvis=Aod_visible(ielem,iline)
          aod_cat=0
          if(aodvis .ge. 30.)then
            aod_cat=1
          elseif(aodvis .ge. 10.)then
            aod_cat=2
          elseif(aodvis .ge. 2.)then
            aod_cat=3
          elseif(aodvis .gt. 0.)then
            aod_cat=4
          endif
          if(aod_cat .ne. 0 )then
            if(modis)then
              visib=.25*visib+.75*(aod_bias(imon,aod_cat) +  &
              aod_slope(imon,aod_cat)*Aod_visible(ielem,iline))
            elseif(abi)then
              visib=.20*visib+.80*(aod_bias(imon,aod_cat) +  &
              aod_slope(imon,aod_cat)*Aod_visible(ielem,iline))
            endif
          endif
!         vibibility in km
!          Visibility(ielem,iline)=min(30.,visib)
          Visibility(ielem,iline)=min(maxvisib,visib) ! new 6/9/2011
          aerosol_visib(ielem,iline)=Visibility(ielem,iline)
          Visib_type(ielem,iline)=0
        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
        if (.not.clear) then
!          now 80 percent threshold
           if(out2(ifog)%MVFR_fog_prob(ielem,iline) > 80. .and.  &
           fog_depth(ielem,iline) .ne. -999.)then
           if(sat%solzen(ielem,iline)<90.and.COT(ielem,iline)>0.0)then
           
            visib_fog=3.0/(COT(ielem,iline)/(fog_depth(ielem,iline)/1000.))

            Fguess_fog_vis(ielem,iline)=visib_fog
            visib_cat=0
            if(visib_fog .ge. 30.)then
              visib_cat=1
            elseif(visib_fog .ge. 10.)then
              visib_cat=2
            elseif(visib_fog .ge. 2.)then
              visib_cat=3
            elseif(visib_fog .gt. 0.)then
              visib_cat=4
            endif
            if(visib_cat .ne. 0 )then
              if(visib_cat .eq. 4)then
                if(modis.or.abi)then
                  visib_fog=.4*visib_fog+.6*(fog_bias(imon,visib_cat) +  &
                  fog_slope(imon,visib_cat)*COT(ielem,iline))
                endif
                if(Visibility(ielem,iline).ne.missing_value_real4)then
                 write(6,*)'should not be in fog',ielem,iline
                endif
!                Visibility(ielem,iline)=min(30.,visib_fog)
                Visibility(ielem,iline)=min(maxvisib,visib_fog) ! 6/9/2011
!                fog_visib(ielem,iline)=min(30.,visib_fog)
                fog_visib(ielem,iline)=min(maxvisib,visib_fog) ! 6/9/2011
                Visib_type(ielem,iline)=1
              else
!               don't expect other categories
                write(6,*)'unexp'
                call flush(6)
                write(6,*)'unexpected fog visibility category at ',ielem, &
                iline,' cat ',visib_cat,' visib_fog ',visib_fog,' COT ', &
                COT(ielem,iline),' fog_depth ',fog_depth(ielem,iline)
!               leave as a missing point for now
              endif
            endif
           endif
          endif
        endif

     end do element_loop
  end do line_loop
  !--- postprocessing
  IF(Sat%Iseg == (Sat%NSeg - 1)) THEN

     DEALLOCATE(COT)
     DEALLOCATE(CWP)
     DEALLOCATE(MOD04_Lon)
     DEALLOCATE(MOD04_Lat)
     DEALLOCATE(MOD06_Lon)
     DEALLOCATE(MOD06_Lat)


  END IF
  deallocate(lwp)
  deallocate(fog_depth)
  write(6,*)'bottom of visib_main INFORMATIONAL'
  call flush(6)


  !--- postprocessing

end SUBROUTINE visib_main
! 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,aod_bias, &
  aod_slope,fog_bias,fog_slope,modis,abi)
implicit none
logical modis,abi
integer(kind=int4) :: nmonbias,ncatbias,nmonbiasi,ncatbiasi
real (kind=real4),dimension(nmonbias,ncatbias) :: aod_bias,aod_slope
real (kind=real4),dimension(nmonbias,ncatbias) :: fog_bias,fog_slope
character *200 filebias,dirbias*(*)
! New for AOD bias and slope
! one for ABI and one for MODIS
! for ABI
if(abi)then
filebias=trim(dirbias)//'/ABI.AOD.LUT.WMO.harmonic.mean.visibility.dat'
endif
! for MODIS
if(modis)then
filebias=trim(dirbias)//'/MODIS.AOD.LUT.WMO.harmonic.mean.visibility.dat'
endif
open(10,file=filebias,form='unformatted')
! read in dimension of bias corrrections
read(10)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)aod_bias,aod_slope
close(10)
if(abi)then
filebias=trim(dirbias)//'/ABI.COT.LUT.WMO.harmonic.mean.visibility.80prob.dat'
elseif(modis)then
filebias=trim(dirbias)//'/MODIS.COT.LUT.WMO.harmonic.mean.visibility.80prob.dat'
endif
!filebias=trim(dirbias)//'/GOES_ASOS_bias_correction.dat'
open(10,file=filebias,form='unformatted')
! read in dimension of bias corrrections
read(10)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)fog_bias,fog_slope
close(10)
return
end SUBROUTINE readvisbiascor
! interpolate from MOD04 resolution to MOD06 resolution( pixel resolution)

! restrictions

! not over a pole or over dateline

subroutine aod04to06(nc6,nr6,lon06,lat06,nc4,nr4,lon04,lat04,nco,nro, &
  aod06,AOD04,ystart,nr6tot)
implicit none
integer(kind=int4) nco,nro,ystart
integer(kind=int4)nc4,nr4,nr6tot
real(kind=real8) lon04(nc4,nr4),lat04(nc4,nr4),AOD04(nc4,nr4)
! need for sat lat lon
!real(kind=real4) lon06(nc6,nr6),lat06(nc6,nr6)
real(kind=real8) lon06(nc6,nr6),lat06(nc6,nr6)
real(kind=real4) aod06(nco,nro)
real(kind=real8) lon04b(nc4,nr4),lat04b(nc4,nr4)
real(kind=real8) dist2,dist2c,drad,pi,cosb
real(kind=real8) lat12,lon12,lat34,lon34,aod12,aod34
real(kind=real8) delta(3),alp0,bet0
real(kind=real8) latbar
integer(kind=int4) ndiff,i1,i2,i3,i4,j1,j2,j3,j4,ipass
! patthern    1          4
!             2          3
! but either 1 or 4 could be north most or 2 and 3 south most
! but 4 is east of 1 and 3 is east or 2 and 1 is north of 2 and 4 is north of 3
!
! set up search pattern for search around nearest box to find box in
parameter (ndiff=3)
integer(kind=int4) nsearch(ndiff)
data nsearch/1,8,16/
integer(kind=int4) idiff(16,ndiff),jdiff(16,ndiff)
data idiff/16*0,-1,0,0,1,-1,-1,1,1,8*0,-2,0,0,0,-2,-2,2,2,-1,-1,1,1,-2,2,2,2/
data jdiff/16*0,0,-1,1,0,-1,1,1,-1,8*0,0,-2,2,2,-1,1,1,-1,-2,2,2,-2,-2,2,2,-2/
data delta/.1,.01,.001/
integer(kind=int4) maxn,i,j,jorder,ii,jj,iii,jjj
integer(kind=int4)iout,nc6,nr6,iorder
integer(kind=int4), allocatable :: num(:,:),ilist(:,:,:),jlist(:,:,:)
integer(kind=int4) nis,nie,njs,nje
integer(kind=int4) maxout,nnum,ib,jb,ibb,jbb
real(kind=real8) alp,bet,latat(0:100,0:100),lonat(0:100,0:100),alpm,betm
real(kind=real8)llat,llon,lonws,lones
parameter (maxout=4) ! maximum boxs can search away from center box
real(kind=real8) dll,a,b,lat,lon
real(kind=real8) lonww(nc4-1,nr4-1),lonee(nc4-1,nr4-1),latnn(nc4-1,nr4-1)
real(kind=real8) latss(nc4-1,nr4-1)
real(kind=real8) lats(nc4-1,nr4-1),latn(nc4-1,nr4-1),lone(nc4-1,nr4-1)
real(kind=real8) lonw(nc4-1,nr4-1)
integer(kind=int4) is,ie,js,je,l,n
integer jbstart
integer entry
save entry
data entry/0/
entry=entry+1
pi=4.*atan(1.)
drad=pi/180.
maxn=1000
dll=0.5
if(lon04(2,1)>lon04(1,1))then
   iorder=1
else
   iorder=-1
endif
if(lat04(1,2)>lat04(1,1))then
   jorder=1
else
   jorder=-1
endif
! set up offset to find pointers to points 1,2,3,4
if(iorder.eq.1)then
  if(jorder.eq.-1)then
    i1=0
    j1=0
    i2=0
    j2=1
    i3=1
    j3=1
    i4=1
    j4=0
  else
    i1=0
    j1=1
    i2=0
    j2=0
    i3=1
    j3=0
    i4=1
    j4=1
  endif
else
  if(jorder.eq.-1)then
    i1=1
    j1=0
    i2=1
    j2=1
    i3=0
    j3=1
    i4=0
    j4=0
  else
    i1=1
    j1=1
    i2=1
    j2=0
    i3=0
    j3=0
    i4=0
    j4=1
  endif
endif
 
if(entry.eq.1)then
 print *,'VISIBILITY aod interpolation (INFORMATION) iorder', &
   iorder,jorder
endif
! final area with non missing lat, lon on MOD04 grid
is=1000000
ie=-100000
js=1000000
je=-1000000
do j=1,nr4
  do i=1,nc4
    if(lat04(i,j) .le. -998. .or. lon04(i,j) .le. -998. )then
!     print *,'missing',i,j
    else
      is=min(i,is)
      js=min(j,js)
      ie=max(i,ie)
      je=max(j,je)
    endif
  end do
end do
! check if a row has zeros
do j=js,je
  do i=is+1,ie-1
    if (lat04(i,j) .eq. 0.0 .and. lon04(i,j) .eq. 0.0)then
      if(lat04(i+1,j) .eq. 0.0 .and. lon04(i+1,j) .eq. 0.0)then
        print *,'cant have row of zeros'
        stop 'zeroes'
      endif
    endif
  end do
end do
nis=100000
njs=100000
nie=-100000
nje=-100000
latnn=1.e32
latss=1.e32
lonee=1.e32
lonww=1.e32
lats=1.e32
latn=1.e32
lone=1.e32
lonw=1.e32
do j=js,je-1
  do i=is,ie-1
!   need a fix if cross +-180 longitude
!   define central point of quadrilateral
    lon04b(i,j)=.25*(lon04(i,j)+lon04(i+1,j)+lon04(i,j+1)+lon04(i+1,j+1))
    lat04b(i,j)=.25*(lat04(i,j)+lat04(i+1,j)+lat04(i,j+1)+lat04(i+1,j+1))
    ii=lon04b(i,j)/dll
    jj=lat04b(i,j)/dll
    
    nis=min(ii,nis)
    njs=min(jj,njs)
    nje=max(jj,nje)
    nie=max(ii,nie)
    ib=i
    jb=j
!   set up out of bounds limits for ib,jb quadrilateral area
    lonww(ib,jb)= &
     min(lon04(i+i1,j+j1),lon04(i+i2,j+j2),lon04(i+i3,j+j3),lon04(i+i4,j+j4))
    lonee(ib,jb)= &
     max(lon04(i+i1,j+j1),lon04(i+i2,j+j2),lon04(i+i3,j+j3),lon04(i+i4,j+j4))
    latnn(ib,jb)= &
     max(lat04(i+i1,j+j1),lat04(i+i2,j+j2),lat04(i+i3,j+j3),lat04(i+i4,j+j4))
    latss(ib,jb)= &
     min(lat04(i+i1,j+j1),lat04(i+i2,j+j2),lat04(i+i3,j+j3),lat04(i+i4,j+j4))
!   set up bounds of inside quasi-rectangular area
    latn(ib,jb)=max(lat04(i+i2,j+j2),lat04(i+i3,j+j3))
    lats(ib,jb)=min(lat04(i+i1,j+j1),lat04(i+i4,j+j4))
    lonw(ib,jb)=min(lon04(i+i3,j+j3),lon04(i+i4,j+j4))
    lone(ib,jb)=max(lon04(i+j1,j+j1),lon04(i+i2,j+j2))
  end do
end do

allocate (num(nis:nie,njs:nje))
num=0
allocate (ilist(nis:nie,njs:nje,maxn),jlist(nis:nie,njs:nje,maxn))
    
do j=js,je-1
  do i=is,ie-1
!   set up hash table i.e. MOD 04 lists to speed up search
    if(lat04(i,j).gt.-998. .and. lon04(i,j).gt.-998.)then
      ii=lon04b(i,j)/dll
      jj=lat04b(i,j)/dll
      if(num(ii,jj)<maxn)then
        num(ii,jj)=num(ii,jj)+1
        ilist(ii,jj,num(ii,jj))=i
        jlist(ii,jj,num(ii,jj))=j
      else
        write(6,*)'num',num(ii,jj),'ii',ii,jj,'lon04b',lon04b(i,j),lat04b(i,j),'i',i,j
        call flush(6)
        print *,'too many'
        stop
      endif
    endif
    
  end do
end do

! now get locations
 do j=1,nr6
  do i=1,nc6
    if(i>5.and.i<nc6-4)then
      ibb=(i-6)/10+1
    else
      ibb=-9999
    endif
    if(j+ystart-1>5.and.j+ystart-1<nr6tot-4)then
      jbb=(j+ystart-1-6)/10+1-max((entry-1)*10-2,0)
    else
      jbb=-9999
    endif
    if(lon06(i,j).gt.-998. .and. lat06(i,j).gt.-998.)then
     nnum=0
     ii=lon06(i,j)/dll
     jj=lat06(i,j)/dll
     ii=max(nis,min(nie,ii))
     jj=max(njs,min(nje,jj))
!    search around
     do iout=1,3
       do n=1,nsearch(iout)
         jjj=jj+jdiff(n,iout)
         iii=ii+idiff(n,iout)
         if(iii<nis.or.iii>nie.or.jjj<njs.or.jjj>nje)go to 1004
         if(num(iii,jjj).ne.0)then
           do l=1,num(iii,jjj)
             jb=jlist(iii,jjj,l)
             ib=ilist(iii,jjj,l)
!            determine if outside ib,jb area
             if(lon06(i,j)<lonww(ib,jb))go to 1021 ! out
             if(lon06(i,j)>lonee(ib,jb))go to 1021 ! out
             if(lat06(i,j)>latnn(ib,jb))goto 1021 ! out
             if(lat06(i,j)<latss(ib,jb))go to 1021 ! out
!            test if insise quasi-rectangular area ib,jb
             if(lat06(i,j)<=lats(ib,jb).and.lat06(i,j)>=latn(ib,jb).and. &
             lon06(i,j)<=lonw(ib,jb).and.lon06(i,j)>=lone(ib,jb))then
               go to 1020 ! inside rectangle
             endif
!            check if with big box
!            check if west of west line i1,i2
             a=lon04(ib+i2,jb+j2)
             b=(lon04(ib+i1,jb+j1)-a)/(lat04(ib+i1,jb+j1)-lat04(ib+i2,jb+j2))
             lon=a+b*(lat06(i,j)-lat04(ib+i2,jb+j2))
             lonws=lon
             if(lon06(i,j)<lon)then
               go to 1021 ! outside
             endif
!            check if south of south line,i2,i3
             a=lat04(ib+i2,jb+j2)
             b=(lat04(ib+i3,jb+j3)-a)/(lon04(ib+i3,jb+j3)-lon04(ib+i2,jb+j2))
             lat=a+b*(lon06(i,j)-lon04(ib+i2,jb+j2))
             if(lat06(i,j)<lat)then
               go to 1021 ! outside
             endif
!            check if east of east line  i3,i4
             a=lon04(ib+i3,jb+j3)
             b=(lon04(ib+i4,jb+j4)-a)/(lat04(ib+i4,jb+j4)-lat04(ib+i3,jb+j3))
             lon=a+b*(lat06(i,j)-lat04(ib+i3,jb+j3))
             lones=lon
             if(lon06(i,j)>lon)then
               go to 1021 ! outside
             endif
!            check if south of north line i1,i4
             a=lat04(ib+i1,jb+j1)
             b=(lat04(ib+i4,jb+j4)-a)/(lon04(ib+i4,jb+j4)-lon04(ib+i1,jb+j1))
             lat=a+b*(lon06(i,j)-lon04(ib+i1,jb+j1))
             if(lat06(i,j)>lat)then
               go to 1021 ! outside
             endif
             go to 1020 ! inside big pixel
   
1021           continue
           end do ! l
         endif
1004     continue
       end do ! n
     end do ! iout
     go to 999
1020  continue
!    box is ib,jb
!    since in box refine search to get alp,bet for 06 point
       dist2=1.e10
       alp0=0.0
       bet0=0.0
       do ipass=1,2
         do jj=0,10
           do ii=0,10
           alp=delta(ipass)*float(ii)+alp0
           bet=delta(ipass)*float(jj)+bet0
           lat12=lat04(ib+i1,jb+j1)*bet+lat04(ib+i2,jb+j2)*(1.-bet)
           lat34=lat04(ib+i4,jb+j4)*bet+lat04(ib+i3,jb+j3)*(1.-bet)
           lon12=lon04(ib+i1,jb+j1)*bet+lon04(ib+i2,jb+j2)*(1.-bet)
           lon34=lon04(ib+i4,jb+j4)*bet+lon04(ib+i3,jb+j3)*(1.-bet)
           latat(ii,jj)=alp*lat34+(1.-alp)*lat12
           lonat(ii,jj)=alp*lon34+(1.-alp)*lon12
           latbar=.5*(lat06(i,j)+latat(ii,jj))*drad
           cosb=cos(latbar)
           dist2c=(lat06(i,j)-latat(ii,jj))**2.+ &
             ((lon06(i,j)-lonat(ii,jj))*cosb)**2.
           if(dist2c<dist2)then
             dist2=dist2c
             alpm=alp
             betm=bet
             llat=latat(ii,jj)
             llon=lonat(ii,jj)
           endif
           end do
         end do
         alp=alpm
         bet=betm
         alp=max(0.0,min(1.0,alp))
         bet=max(0.0,min(1.0,bet))
         alp0=alp-0.5*delta(ipass)
         bet0=bet-0.5*delta(ipass)
         alp0=max(0.0,min(1.0,alp0))
         bet0=max(0.0,min(1.0,bet0))
       end do
       
       if(AOD04(ib+i1,jb+j1)>=0.0.and.AOD04(ib+i2,jb+j2)>=0.0.and. &
          AOD04(ib+i3,jb+j3)>=0.0.and.AOD04(ib+i4,jb+j4)>=0.0)then
         aod12=AOD04(ib+i1,jb+j1)*bet+AOD04(ib+i2,jb+j2)*(1.-bet)
         aod34=AOD04(ib+i4,jb+j4)*bet+AOD04(ib+i3,jb+j3)*(1.-bet)
         aod06(i,j)=alp*aod34+(1.-alp)*aod12

       endif
999  continue
   endif
  end do ! i
 end do ! j
deallocate (num,ilist,jlist)
return
end subroutine aod04to06

end MODULE visibmod 
