source: trunk/LMDZ.TITAN/libf/phytitan/disr_haze.F90 @ 1648

subroutine disr_haze(dz,press,wno,taeros,ssa,cbar)
use datafile_mod, only: datadir
IMPLICIT NONE

real*8,intent(in)   :: dz,press,wno
real*8,intent(inout):: taeros,ssa,cbar

!---------------------------
! Attention !!
! taeros is the integrated extinction over the layer
! (extinction * thickness of layer)
!---------------------------

integer             :: i,j,iw,ip,ierr
real*8              :: wln,factw,factp
integer,parameter   :: nbwl_PL=328,nblev_PL=162
real*8,save         :: ext_PL(nblev_PL,nbwl_PL),ssa_PL(nblev_PL,nbwl_PL)
real*8,save         :: asf_PL(nblev_PL,nbwl_PL)
real*8,save         :: wl_PL(nbwl_PL),press_PL(nblev_PL)
logical,save        :: firstcall=.true.
character(len=15)   :: dummy

if (firstcall) then
   print*,"DISR HAZE - PANAYOTIS LAVVAS"


! read PL table   
! wl_PL in nm
! press_PL in Pa
   open(11,file=TRIM(datadir)//'/hazetable_PL_original.dat',status="old",iostat=ierr)
   read(11,*) dummy,wl_PL
   do i=1,nblev_PL
     read(11,*) press_PL(i),ext_PL(i,:)  ! in cm-1
   enddo
   do i=1,nblev_PL
     read(11,*) press_PL(i),ssa_PL(i,:)
   enddo
   do i=1,nblev_PL
     read(11,*) press_PL(i),asf_PL(i,:)
   enddo
   close(11)
   ! convert press_PL into millibar for comparison to press in the generic (wtf press in mbar but ok ... to be modified at some point anyway)
   press_PL(:)=press_PL(:)*1E-2

   firstcall=.false.
endif

! convert wno (in cm-1) into wln (nm)
wln=1E7/wno

! interpolate the needed values from the table

iw=1
do i=2,nbwl_PL
  if(wln.gt.wl_PL(i)) then
    iw=i
  endif
enddo

ip=1
do j=2,nblev_PL
  if(press.lt.press_PL(j)) then
    ip=j
  endif
enddo

if(iw.ne.nbwl_PL) then
  factw=(wln-wl_PL(iw))/(wl_PL(iw+1)-wl_PL(iw))
endif
if(ip.ne.nblev_PL) then
  factp=(log10(press)-log10(press_PL(ip)))/(log10(press_PL(ip+1))-log10(press_PL(ip)))
else
!  factp=(log10(press)-log10(press_PL(nblev_PL)))/(log10(1.e-10)-log10(press_PL(nblev_PL)))
  factp=0.
endif

if ((iw.ne.nbwl_PL).and.(ip.ne.nblev_PL)) then
  taeros= ext_PL(ip,iw)  *(1-factw)*(1-factp)+ext_PL(ip+1,iw)  *(1-factw)*factp &
         +ext_PL(ip,iw+1)*factw    *(1-factp)+ext_PL(ip+1,iw+1)*factw    *factp
  ssa   = ssa_PL(ip,iw)  *(1-factw)*(1-factp)+ssa_PL(ip+1,iw)  *(1-factw)*factp &
         +ssa_PL(ip,iw+1)*factw    *(1-factp)+ssa_PL(ip+1,iw+1)*factw    *factp
  cbar  = asf_PL(ip,iw)  *(1-factw)*(1-factp)+asf_PL(ip+1,iw)  *(1-factw)*factp &
         +asf_PL(ip,iw+1)*factw    *(1-factp)+asf_PL(ip+1,iw+1)*factw    *factp
else if ((iw.eq.nbwl_PL).and.(ip.ne.nblev_PL)) then
  taeros= ext_PL(ip,iw)*(1-factp)+ext_PL(ip+1,iw)*factp 
  ssa   = ssa_PL(ip,iw)*(1-factp)+ssa_PL(ip+1,iw)*factp 
  cbar  = asf_PL(ip,iw)*(1-factp)+asf_PL(ip+1,iw)*factp
else if ((iw.ne.nbwl_PL).and.(ip.eq.nblev_PL)) then
  taeros= ext_PL(ip,iw)  *(1-factw)*(1-factp) &
         +ext_PL(ip,iw+1)*factw    *(1-factp)
  ssa   = ssa_PL(ip,iw)  *(1-factw)*(1-factp) &
         +ssa_PL(ip,iw+1)*factw    *(1-factp)
  cbar  = asf_PL(ip,iw)  *(1-factw)*(1-factp) &
         +asf_PL(ip,iw+1)*factw    *(1-factp)
else if ((iw.eq.nbwl_PL).and.(ip.eq.nblev_PL)) then
  taeros= ext_PL(ip,iw)*(1-factp) 
  ssa   = ssa_PL(ip,iw)*(1-factp) 
  cbar  = asf_PL(ip,iw)*(1-factp)
endif

taeros=taeros*dz*1.E2 ! ext in cm-1 * thickness in m * 1E2

end subroutine disr_haze