source: trunk/LMDZ.GENERIC/libf/phystd/interpolateH2Ocont.F90 @ 603

     subroutine interpolateH2Ocont(wn,temp,presS,presF,abcoef,firstcall)

!==================================================================
!     
!     Purpose
!     -------
!     Calculates the H2O continuum opacity, using a lookup table from
!     Clough (2005)
!
!     Authors
!     -------
!     R. Wordsworth (2011)
!     
!==================================================================

      use datafile_mod, only: datadir
      implicit none

      ! input
      double precision wn                 ! wavenumber             (cm^-1)
      double precision temp               ! temperature            (Kelvin)
      double precision presS              ! self-pressure          (Pascals)
      double precision presF              ! foreign (air) pressure (Pascals)

      ! output
      double precision abcoef             ! absorption coefficient (m^-1)

      integer nS,nT
      parameter(nS=1001)
      parameter(nT=11)

      double precision kB
      parameter(kB=1.3806488e-23)

      double precision amagatS, amagatF, abcoefS, abcoefF, Nmolec
      double precision wn_arr(nS)
      double precision temp_arr(nT)
      double precision abs_arrS(nS,nT)
      double precision abs_arrF(nS,nT)
      double precision data_tmp(nT)

      integer k
      logical firstcall

      save wn_arr, temp_arr, abs_arrS, abs_arrF

      character*100 dt_file
      integer strlen,ios

      amagatS=(273.15/temp)*(presS/101325.0)
      amagatF=(273.15/temp)*(presF/101325.0)

      if(firstcall)then ! called by sugas_corrk only
         print*,'----------------------------------------------------'
         print*,'Initialising H2O continuum from MT_CKD data...'

!     1.1 Open the ASCII files

         ! nu array
         dt_file=TRIM(datadir)//'/continuum_data/H2O_CONT_NU.dat'
         open(33,file=dt_file,form='formatted',status='old',iostat=ios)
         if (ios.ne.0) then        ! file not found
           write(*,*) 'Error from interpolateH2O_cont'
           write(*,*) 'Data file ',trim(dt_file),' not found.'
           write(*,*)'Check that your path to datagcm:',trim(datadir)
           write(*,*)' is correct. You can change it in callphys.def with:'
           write(*,*)' datadir = /absolute/path/to/datagcm'
           write(*,*)' Also check that there is a continuum_data/H2O_CONT_NU.dat there.'
           call abort
         else
            do k=1,nS
               read(33,*) wn_arr(k)
            enddo
         endif
         close(33)

         ! self broadening
         dt_file=TRIM(datadir)//'/continuum_data/H2O_CONT_SELF.dat'
         open(34,file=dt_file,form='formatted',status='old',iostat=ios)
         if (ios.ne.0) then        ! file not found
           write(*,*) 'Error from interpolateH2O_cont'
           write(*,*) 'Data file ',trim(dt_file),' not found.'
           write(*,*)'Check that your path to datagcm:',trim(datadir)
           write(*,*)' is correct. You can change it in callphys.def with:'
           write(*,*)' datadir = /absolute/path/to/datagcm'
           write(*,*)' Also check that there is a continuum_data/H2O_CONT_SELF.dat there.'
           call abort
         else
            do k=1,nS
               read(34,*) data_tmp
               abs_arrS(k,1:nT)=data_tmp(1:nT)
            end do
         endif
         close(34)

         ! foreign (N2+O2+Ar) broadening
         dt_file=TRIM(datadir)//'/continuum_data/H2O_CONT_FOREIGN.dat'
         open(35,file=dt_file,form='formatted',status='old',iostat=ios)
         if (ios.ne.0) then        ! file not found
           write(*,*) 'Error from interpolateH2O_cont'
           write(*,*) 'Data file ',trim(dt_file),' not found.'
           write(*,*)'Check that your path to datagcm:',trim(datadir)
           write(*,*)' is correct. You can change it in callphys.def with:'
           write(*,*)' datadir = /absolute/path/to/datagcm'
           write(*,*)' Also check that there is a continuum_data/H2O_CONT_FOREIGN.dat there.'
           call abort
         else
            do k=1,nS
               read(35,*) data_tmp
               abs_arrF(k,1:nT)=data_tmp(1:nT)
            end do
         endif
         close(35)

         temp_arr(1)  = 200.
         temp_arr(2)  = 250.
         temp_arr(3)  = 300.
         temp_arr(4)  = 350.
         temp_arr(5)  = 400.
         temp_arr(6)  = 450.
         temp_arr(7)  = 500.
         temp_arr(8)  = 550.
         temp_arr(9)  = 600.
         temp_arr(10) = 650.
         temp_arr(11) = 700.

         print*,'interpolateH2Ocont: At wavenumber ',wn,' cm^-1'
         print*,'   temperature ',temp,' K'
         print*,'   H2O pressure ',presS,' Pa'
         print*,'   air pressure ',presF,' Pa'

         call bilinearH2Ocont(wn_arr,temp_arr,abs_arrS,wn,temp,abcoefS)
         print*,'the self absorption is ',abcoefS,' cm^2 molecule^-1'

         call bilinearH2Ocont(wn_arr,temp_arr,abs_arrF,wn,temp,abcoefF)
         print*,'the foreign absorption is ',abcoefF,' cm^2 molecule^-1'

         print*,'We have ',amagatS,' amagats of H2O vapour'
         print*,'and ',amagatF,' amagats of air'

         abcoef = abcoefS*amagatS + abcoefF*amagatF ! Eq. (15) in Clough (1989)

         Nmolec = (presS+presF)/(kB*temp)

         print*,'Total number of molecules per m^3 is',Nmolec

         abcoef = abcoef*Nmolec/(100.0**2) ! convert to m^-1
         abcoef = abcoef*(presS/(presF+presS))     ! take H2O mixing ratio into account

         print*,'So the total absorption is ',abcoef,' m^-1'
         print*,'And optical depth / km : ',1000.0*abcoef

      else

         call bilinearH2Ocont(wn_arr,temp_arr,abs_arrS,wn,temp,abcoefS)
         call bilinearH2Ocont(wn_arr,temp_arr,abs_arrF,wn,temp,abcoefF)

         abcoef = abcoefS*amagatS + abcoefF*amagatF ! Eq. (15) in Clough (1989)

         !abcoef = (presS/100000.0)*3.0e-24
         !print*,'Matsui TEST'

         Nmolec = (presS+presF)/(kB*temp)  ! sure this is correct??
         abcoef = abcoef*Nmolec/(100.0**2) ! convert to m^-1
         abcoef = abcoef*(presS/(presF+presS))     ! take H2O mixing ratio into account


         ! unlike for Rayleigh scattering, we do not currently weight by the BB function
         ! however our bands are normally thin, so this is no big deal.


      endif

      return
    end subroutine interpolateH2Ocont




!-------------------------------------------------------------------------
      subroutine bilinearH2Ocont(x_arr,y_arr,f2d_arr,x_in,y_in,f)
!     Necessary for interpolation of continuum data

      implicit none

      integer nX,nY,i,j,a,b
      parameter(nX=1001)
      parameter(nY=11)

      real*8 x_in,y_in,x,y,x1,x2,y1,y2
      real*8 f,f11,f12,f21,f22,fA,fB
      real*8 x_arr(nX)
      real*8 y_arr(nY)
      real*8 f2d_arr(nX,nY)
      
      integer strlen
      character*100 label
      label='subroutine bilinear'

      x=x_in
      y=y_in

!     1st check we're within the wavenumber range
      if ((x.lt.x_arr(2)).or.(x.gt.x_arr(nX-2))) then
         f=0.0D+0
         return
      else
         
!     in the x (wavenumber) direction 1st
         i=1
 10      if (i.lt.(nX+1)) then
            if (x_arr(i).gt.x) then
               x1=x_arr(i-1)
               x2=x_arr(i)
               a=i-1
               i=9999
            endif
            i=i+1
            goto 10
         endif
      endif
      
      if ((y.lt.y_arr(1)).or.(y.gt.y_arr(nY))) then
         write(*,*) 'Warning from bilinearH2Ocont:'
         write(*,*) 'Outside continuum temperature range!'
         if(y.lt.y_arr(1))then
            y=y_arr(1)+0.01
         endif
         if(y.gt.y_arr(nY))then
            y=y_arr(nY)-0.01
         endif
      else

!     in the y (temperature) direction 2nd
         j=1
 20      if (j.lt.(nY+1)) then
            if (y_arr(j).gt.y) then
               y1=y_arr(j-1)
               y2=y_arr(j)
               b=j-1
               j=9999
            endif
            j=j+1
            goto 20
         endif
      endif
      
      f11=f2d_arr(a,b)
      f21=f2d_arr(a+1,b)
      f12=f2d_arr(a,b+1)
      f22=f2d_arr(a+1,b+1)
      
!     1st in x-direction
      fA=f11*(x2-x)/(x2-x1)+f21*(x-x1)/(x2-x1)
      fB=f12*(x2-x)/(x2-x1)+f22*(x-x1)/(x2-x1)
      
!     then in y-direction
      f=fA*(y2-y)/(y2-y1)+fB*(y-y1)/(y2-y1)

      return
    end subroutine bilinearH2Ocont