text:
      program atmos_gaussian
      implicit none
c
c --- Usage: echo JH | atmos_gaussian
c ---   print the 2*JH latitudes of the specified gaussian grid
c
c --- A "Gaussian" grid, is regular in longitude and
c ---   almost regular in latitude (Hortal and Simmons, 1991).
c --- https://en.wikipedia.org/wiki/Gaussian_grid
c
      integer              :: jh,j
      real*8               :: pscale
      real*4,  allocatable :: plat(:)
      real*8,  allocatable :: slat(:)
c
      read(5,*) jh
c
      allocate( plat(2*jh) )
      allocate( slat(  jh) )
c
      call glat8(jh,slat)
      pscale = 180.d0/acos(-1.d0)
      do j= 1,jh
        plat(2*jh+1-j) =  pscale*asin(slat(j))
        plat(       j) = -plat(2*jh+1-j)
      enddo !j
      do j= 1,2*jh
        write(6,'(f10.5)') plat(j)
      enddo
      end
      SUBROUTINE GLAT8(JH,SLAT)
      IMPLICIT NONE
C
      INTEGER JH
      REAL*8  SLAT(JH)
C
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:    GLAT        COMPUTE GAUSSIAN LATITUDE FUNCTIONS
C   PRGMMR: IREDELL          ORG: W/NMC23     DATE: 92-10-31
C
C ABSTRACT: COMPUTES SINES OF GAUSSIAN LATITUDE BY ITERATION.
C
C PROGRAM HISTORY LOG:
C   91-10-31  MARK IREDELL
C
C USAGE:    CALL GLAT(JH,SLAT)
C
C   INPUT ARGUMENT LIST:
C     JH       - INTEGER NUMBER OF GAUSSIAN LATITUDES IN A HEMISPHERE
C
C   OUTPUT ARGUMENT LIST:
C     SLAT     - REAL (JH) SINES OF (POSITIVE) GAUSSIAN LATITUDE
C
C ATTRIBUTES:
C   LANGUAGE: CRAY FORTRAN
C
C$$$
C
      INTEGER   JBZ
      PARAMETER(JBZ=50)
      REAL*8    PI,C,EPS
      PARAMETER(PI=3.14159265358979,C=(1.-(2./PI)**2)*0.25,EPS=1.E-14)
C
      INTEGER ITS,J,N
      REAL*8  PKM2,R,SP,SPMAX
      REAL*8  PK(JH),PKM1(JH)
C
      REAL*8  BZ(JBZ)
      DATA BZ        / 2.4048255577,  5.5200781103,
     $  8.6537279129, 11.7915344391, 14.9309177086, 18.0710639679,
     $ 21.2116366299, 24.3524715308, 27.4934791320, 30.6346064684,
     $ 33.7758202136, 36.9170983537, 40.0584257646, 43.1997917132,
     $ 46.3411883717, 49.4826098974, 52.6240518411, 55.7655107550,
     $ 58.9069839261, 62.0484691902, 65.1899648002, 68.3314693299,
     $ 71.4729816036, 74.6145006437, 77.7560256304, 80.8975558711,
     $ 84.0390907769, 87.1806298436, 90.3221726372, 93.4637187819,
     $ 96.6052679510, 99.7468198587, 102.888374254, 106.029930916,
     $ 109.171489649, 112.313050280, 115.454612653, 118.596176630,
     $ 121.737742088, 124.879308913, 128.020877005, 131.162446275,
     $ 134.304016638, 137.445588020, 140.587160352, 143.728733573,
     $ 146.870307625, 150.011882457, 153.153458019, 156.295034268 /
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  ESTIMATE LATITUDES USING BESSEL FUNCTION
      R=1./SQRT((2*JH+0.5)**2+C)
      DO J=1,MIN(JH,JBZ)
        SLAT(J)=COS(BZ(J)*R)
      ENDDO
      DO J=JBZ+1,JH
        SLAT(J)=COS((BZ(JBZ)+(J-JBZ)*PI)*R)
      ENDDO
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  CONVERGE UNTIL ALL SINES OF GAUSSIAN LATITUDE ARE WITHIN EPS
      DO 110 ITS= 1,999
        SPMAX=0.
        DO J=1,JH
          PKM1(J)=1.
          PK(J)=SLAT(J)
        ENDDO
        DO N=2,2*JH
          DO J=1,JH
            PKM2=PKM1(J)
            PKM1(J)=PK(J)
            PK(J)=((2*N-1)*SLAT(J)*PKM1(J)-(N-1)*PKM2)/N
          ENDDO
        ENDDO
        DO J=1,JH
          SP=PK(J)*(1.-SLAT(J)**2)/(2*JH*(PKM1(J)-SLAT(J)*PK(J)))
          SLAT(J)=SLAT(J)-SP
          SPMAX=MAX(SPMAX,ABS(SP))
        ENDDO
        IF     (SPMAX.LE.EPS) THEN
          GOTO 1110
        ENDIF
  110 CONTINUE
 1110 CONTINUE
      RETURN
C     END OF GLAT8.
      END

ext: f

max_stars_repo_name: TillRasmussen/HYCOM-tools

max_issues_repo_name: TillRasmussen/HYCOM-tools

max_forks_repo_name: TillRasmussen/HYCOM-tools


SUBROUTINE rmax ( a, n, real_max, imax )
!
!  Purpose:
!    To find the maximum value in an array, and the location
!    of that value in the array.
!
IMPLICIT NONE

! Data dictionary: declare calling parameter types & definitions
INTEGER, INTENT(IN) :: n             ! No. of vals in array a.
REAL, INTENT(IN), DIMENSION(n) :: a  ! Input data.
REAL, INTENT(OUT) :: real_max        ! Maximum value in a.
INTEGER, INTENT(OUT) :: imax         ! Location of max value.

! Data dictionary: declare local variable types & definitions
INTEGER :: i                         ! Index variable

! Initialize the maximum value to first value in array.
real_max = a(1)
imax = 1

! Find the maximum value.
DO i = 2, n
   IF ( a(i) > real_max ) THEN
      real_max = a(i)
      imax = i
   END IF
END DO

END SUBROUTINE rmax

SUBROUTINE rmin ( a, n, real_min, imin )
!
!  Purpose:
!    To find the minimum value in an array, and the location
!    of that value in the array.
!
IMPLICIT NONE

! Data dictionary: declare calling parameter types & definitions
INTEGER, INTENT(IN) :: n             ! No. of vals in array a.
REAL, INTENT(IN), DIMENSION(n) :: a  ! Input data.
REAL, INTENT(OUT) :: real_min        ! Minimum value in a.
INTEGER, INTENT(OUT) :: imin         ! Location of min value.

! Data dictionary: declare local variable types & definitions
INTEGER :: i                         ! Index variable

! Initialize the minimum value to first value in array.
real_min = a(1)
imin = 1

! Find the minimum value.
DO I = 2, n
   IF ( a(i) < real_min ) THEN
      real_min  = a(i)
      imin = i
   END IF
END DO

END SUBROUTINE rmin

SUBROUTINE ave_sd ( a, n, ave, std_dev, error )
!
!  Purpose:
!    To calculate the average and standard deviation of an array.
!
IMPLICIT NONE

! Data dictionary: declare calling parameter types & definitions
INTEGER, INTENT(IN) :: n             ! No. of vals in array a.
REAL, INTENT(IN), DIMENSION(n) :: a  ! Input data.
REAL, INTENT(OUT) :: ave             ! Average of a.
REAL, INTENT(OUT) :: std_dev         ! Standard deviation.
INTEGER, INTENT(OUT) :: error        ! Flag: 0 -- no error
                                     !       1 -- sd invalid
                                     !       2 -- ave & sd invalid

! Data dictionary: declare local variable types & definitions
INTEGER :: i                         ! Loop index
REAL :: sum_x                        ! Sum of input values
REAL :: sum_x2                       ! Sum of input values squared

! Initialize the sums to zero.
sum_x  = 0.
sum_x2 = 0.

! Accumulate sums.
DO I = 1, n
   sum_x  = sum_x + a(i)
   sum_x2 = sum_x2 + a(i)**2
END DO

! Check to see if we have enough input data.
IF ( n >= 2 ) THEN ! we have enough data

   ! Calculate the mean and standard deviation
   ave   = sum_x / REAL(n)
   std_dev = SQRT( (REAL(n) * sum_x2 - sum_x**2) &
           / (REAL(n) * REAL(n-1)) )
   error = 0

ELSE IF ( n == 1 ) THEN ! no valid std_dev

   ave   = sum_x
   std_dev = 0.              ! std_dev invalid
   error = 1

ELSE

   ave = 0.                  ! ave invalid
   std_dev = 0.              ! std_dev invalid
   error = 2

END IF
END SUBROUTINE ave_sd

SUBROUTINE median ( a, n, med )
!
!  Purpose:
!    To calculate the median value of an array.
!
IMPLICIT NONE

! Data dictionary: declare calling parameter types & definitions
INTEGER, INTENT(IN) :: n             ! No. of vals in array a.
REAL, INTENT(IN), DIMENSION(n) :: a  ! Input data.
REAL, INTENT(OUT) :: med             ! Median value of a.

! Sort the data into ascending order.
CALL sort ( a, n )

! Get median.
IF ( MOD(n,2) == 0 ) THEN
   med = ( a(n/2) + a(n/2+1) ) / 2.
ELSE
   med = a(n/2+1)
END IF
END SUBROUTINE median


! written by jxzou at 20210328: adjust the orders of d,f,g, etc functions in
!  .fch(k) file, to the orders in QM4D

! Note:
! 1) QM4D supports only Cartesian-type basis set (6D 10F)
! 2) For ECP integrals, QM4D use spherical harmonic functions. For example,
!    energies of RHF/LANL2DZ for Zn atom is different in QM4D vs Gaussian,
!    but MOs are the same.

! Restrictions:
! 1) Atoms belong to an element must have only one basis set
!    C1 cc-pVTZ with C2 cc-pVDZ is not supported
! 2) Background charges are not stored in .fch(k) file, so the generated .inp
!    file will not contain background charges, neither
! 3) Currently only HF/DFT wavefunction is supported

! The 'Shell types' array in Gaussian .fch file:
!
!   Spherical     |     Cartesian
! -5,-4,-3,-2,-1, 0, 1, 2, 3, 4, 5
!  H  G  F  D  L  S  P  D  F  G  H
!
! 'L' is 'SP' in Pople-type basis sets

program main
 implicit none
 integer :: i
 integer, parameter :: iout = 6
 character(len=4) :: str
 character(len=240) :: fname
 logical :: binary

 i = iargc()
 if(i<1 .or. i>2) then
  write(iout,'(/,1X,A)') 'ERROR in subroutine fch2qm4d: wrong command line arguments!'
  write(iout,'(1X,A)') 'Example 1 (for HF/DFT): fch2qm4d a.fch'
  write(iout,'(1X,A)') 'Example 2 (for HF/DFT): fch2qm4d a.fch -xml'
  write(iout,'(1X,A,/)') 'Example 3 (for HF/DFT): fch2qm4d a.fch -bin'
  stop
 end if

 fname = ' '
 call getarg(1, fname)
 call require_file_exist(fname)

 binary = .false.
 ! False: use .xml ASCII text file
 ! True : use .bin binary file

 if(i == 2) then
  str = ' '
  call getarg(2, str)
  select case(str)
  case('-xml')
  case('-bin')
   binary = .true.
  case default
   write(iout,'(/,A)') 'ERROR in subroutine fch2qm4d: the 2nd command line&
                      & argument is wrong.'
   write(iout,'(A)') "Only '-xml' or '-bin' is accepted. But got '"//str//"'"
   stop
  end select
 end if

 call fch2qm4d(fname, binary)
 stop
end program main

! read the MOs in .fch(k) file and adjust its d,f,g, etc. functions order
!  of Gaussian to that of QM4D
subroutine fch2qm4d(fchname, binary)
 use fch_content, only: check_uhf_in_fch, read_fch, nif, nbf, ncontr, shell_type,&
  natom, coor, elem, na, nb, alpha_coeff, beta_coeff
 implicit none
 integer :: i, j, k, fid
 integer :: n6dmark, n10fmark
 integer, parameter :: iout = 6
 integer, allocatable :: d_mark(:), f_mark(:)
 character(len=240) :: inpname ! QM4D input file
 character(len=240) :: xyzname ! .xyz file to hold Cartesian coordinates
 character(len=240) :: xmlname ! .xml file to hold density matrix
 character(len=240) :: binname ! .bin file to hold MOs (binary)
 character(len=240), intent(in) :: fchname
 real(kind=8), allocatable :: coeff(:,:), alpha_occ(:,:), beta_occ(:,:)
 real(kind=8), allocatable :: alpha_dm(:,:), beta_dm(:,:)
 logical :: uhf
 logical, intent(in) :: binary

 if(binary) then
  write(iout,'(A)') 'ERROR in subroutine fch2qm4d: currently binary=.True.&
                   & is not supported.'
  write(iout,'(A)') "Please use '-xml'."
  stop
 end if

 call check_uhf_in_fch(fchname, uhf) ! determine whether UHF
 call read_fch(fchname, uhf)         ! read content in .fch

 ! check if any Cartesian functions
 if( ANY(shell_type < -1) ) then
  write(iout,'(A)') 'ERROR in subroutine fch2qm4d: spherical harmonic functions&
                   & detected in file '//TRIM(fchname)//'.'
  write(iout,'(A)') "QM4D currently supports only Cartesian functions. You need to add&
                  & '6D 10F' keywords in Gaussian."
  stop
 end if

 ! check any angular momentum higher than F
 if( ANY(shell_type > 3) ) then
  write(iout,'(A)') 'ERROR in subroutine fch2qm4d: it seems that QM4D does not&
                   & support G, H functions.'
  write(iout,'(A)') 'You can use a smaller basis set (e.g. cc-pVQZ -> cc-pVTZ).'
  stop
 end if

 i = index(fchname,'.fch',back=.true.)
 inpname = fchname(1:i-1)//'.inp'
 xyzname = fchname(1:i-1)//'.xyz'
 xmlname = fchname(1:i-1)//'_0.xml'
 binname = fchname(1:i-1)//'_0.bin'

 ! create the .xyz file
 call write_xyzfile(natom, coor, elem, xyzname)
 deallocate(coor)

 ! create the QM4D input file
 call write_qm4d_inp(inpname, uhf, binary)

 ! create basis set files *.gen for each element
 call write_basis_files_for_each_elem()

 if(uhf) then
  allocate(coeff(nbf,2*nif), source=0d0)
  coeff(1:nbf,1:nif) = alpha_coeff
  coeff(1:nbf,nif+1:2*nif) = beta_coeff
  nif = 2*nif
 else
  allocate(coeff(nbf,nif), source=alpha_coeff)
 end if
 deallocate(alpha_coeff)

! then we adjust the basis functions in each MO according to the type of basis functions
 n6dmark = 0
 n10fmark = 0
 allocate(d_mark(ncontr), f_mark(ncontr))
 d_mark = 0
 f_mark = 0
 nbf = 0
 do i = 1, ncontr, 1
  select case(shell_type(i))
  case( 0)   ! S
   nbf = nbf + 1
  case( 1)   ! 3P
   nbf = nbf + 3
  case(-1)   ! SP or L
   nbf = nbf + 4
  case( 2)   ! 6D
   n6dmark = n6dmark + 1
   d_mark(n6dmark) = nbf + 1
   nbf = nbf + 6
  case( 3)   ! 10F
   n10fmark = n10fmark + 1
   f_mark(n10fmark) = nbf + 1
   nbf = nbf + 10
  case default
   write(iout,'(A)') 'ERROR in subroutine fch2qm4d: shell_type out of range!'
   write(iout,'(2(A,I0))') 'i=', i, ', shell_type(i)=', shell_type(i)
   stop
  end select
 end do ! for i

 deallocate(shell_type)

 ! adjust the order of d, f, etc. functions
 do i = 1, n6dmark, 1
  call fch2qm4d_permute_6d(nif,coeff(d_mark(i):d_mark(i)+5,:))
 end do
 do i = 1, n10fmark, 1
  call fch2qm4d_permute_10f(nif,coeff(f_mark(i):f_mark(i)+9,:))
 end do
! adjustment finished

 deallocate(d_mark, f_mark)

 ! initilize occupation number arrays
 if(uhf) then
  nif = nif/2
  allocate(alpha_occ(nif,nif), source=0d0)
  forall(i=1:na) alpha_occ(i,i) = 1d0
  allocate(beta_occ(nif,nif), source=0d0)
  forall(i=1:nb) beta_occ(i,i) = 1d0
 else
  allocate(alpha_occ(nif,nif), source=0d0)
  forall(i=1:nb) alpha_occ(i,i) = 2d0
  if(na > nb) forall(i=nb+1:na) alpha_occ(i,i) = 1d0
 end if

 ! calculate alpha density (and beta density, if any)
 alpha_coeff = coeff(:,1:nif)
 allocate(alpha_dm(nbf,nbf))
 alpha_dm = MATMUL(MATMUL(alpha_coeff,alpha_occ),TRANSPOSE(alpha_coeff))
 deallocate(alpha_coeff)
 if(uhf) then
  beta_coeff = coeff(:,nif+1:2*nif)
  allocate(beta_dm(nbf,nbf))
  beta_dm = MATMUL(MATMUL(beta_coeff,beta_occ),TRANSPOSE(beta_coeff))
  deallocate(beta_coeff)
 end if
 deallocate(coeff)

 ! print alpha density and beta density into _0.xml file
 open(newunit=fid,file=TRIM(xmlname),status='replace')
 write(fid,'(A)') "<?xml version=""1.0"" encoding=""ISO-8859-1""?>"
 write(fid,'(A)') '<QM4D>'
 write(fid,'(A)') ' <!--The symmetric matrix is saved as the lower triangle format.-->'
 write(fid,'(2(A,I6),A)') " <Alpha_Density_Matrix Dimension=""", nbf, """ Data=""",&
                          nbf*(nbf+1)/2,""">"
 k = 0
 do i = 1, nbf, 1
  do j = 1, i, 1
   k = k + 1
   write(fid,'(A,I6,A,1X,ES23.16,A)') "  <Data id=""",k,""">",alpha_dm(j,i),'</Data>'
  end do ! for j
 end do ! for i

 deallocate(alpha_dm)
 write(fid,'(A)') ' </Alpha_Density_Matrix>'

 write(fid,'(A)') ' <!--The full matrix is saved as nRow*nCol.-->'
 write(fid,'(2(A,I6),A)') " <Alpha_Occup_No Row=""     1"" Col=""",nif,""" Data=""",nif,""">"
 do i = 1, nif, 1
  write(fid,'(A,I6,A,1X,ES23.16,A)') "  <Data id=""",i,""">",alpha_occ(i,i),'</Data>'
 end do ! for i
 deallocate(alpha_occ)
 write(fid,'(A)') '  </Alpha_Occup_No>'

 if(uhf) then
  write(fid,'(A)') ' <!--The symmetric matrix is saved as the lower triangle format.-->'
  write(fid,'(2(A,I6),A)') " <Beta_Density_Matrix Dimension=""", nbf, """ Data=""",&
                           nbf*(nbf+1)/2,""">"
  k = 0
  do i = 1, nbf, 1
   do j = 1, i, 1
    k = k + 1
    write(fid,'(A,I6,A,1X,ES23.16,A)') "  <Data id=""",k,""">",beta_dm(j,i),'</Data>'
   end do ! for j
  end do ! for i
  deallocate(beta_dm)
  write(fid,'(A)') ' </Beta_Density_Matrix>'

  write(fid,'(A)') ' <!--The full matrix is saved as nRow*nCol.-->'
  write(fid,'(2(A,I6),A)') " <Beta_Occup_No Row=""     1"" Col=""",nif,""" Data=""",nif,""">"
  do i = 1, nif, 1
   write(fid,'(A,I6,A,1X,ES23.16,A)') "  <Data id=""",i,""">",beta_occ(i,i),'</Data>'
  end do ! for i
  deallocate(beta_occ)
  write(fid,'(A)') '  </Beta_Occup_No>'
 end if

 write(fid,'(A)') '</QM4D>'
 close(fid)
 return
end subroutine fch2qm4d

! write/create a QM4D .inp file
subroutine write_qm4d_inp(inpname, uhf, binary)
 use fch_content, only: charge, mult, natom, elem, LPSkip
 implicit none
 integer :: i, fid
 character(len=240) :: xml0name, xmlname, binname, xyzname
 character(len=240), intent(in) :: inpname
 logical, intent(in) :: uhf, binary

 i = index(inpname, '.inp', back=.true.)
 xml0name = inpname(1:i-1)//'_0.xml'
 xmlname = inpname(1:i-1)//'.xml'
 binname = inpname(1:i-1)//'.bin'
 xyzname = inpname(1:i-1)//'.xyz'

 open(newunit=fid,file=TRIM(inpname),status='replace')
 write(fid,'(A)') '$qm'
 write(fid,'(A)') 'xyz    '//TRIM(xyzname)
 if(uhf) then
  write(fid,'(A)') 'spin   2'
 else
  write(fid,'(A)') 'spin   1'
 end if
 write(fid,'(A,I0)') 'charge ', charge
 write(fid,'(A,I0)') 'mult   ', mult
 write(fid,'(A)') 'iter   128'
 write(fid,'(A)') 'method hf'
 write(fid,'(A)') 'etol   1.e-08'
 if(binary) then
  write(fid,'(A)') 'guess  read '//TRIM(binname)
 else
  write(fid,'(A)') 'guess  read '//TRIM(xml0name)
 end if
 write(fid,'(A)') 'diis   6 1000.0'

 do i = 1, natom, 1
  if(i > 1) then
   if(ANY(elem(1:i-1) == elem(i))) cycle
  end if

  if(allocated(LPSkip)) then
   if(LPSkip(i) == 0) then
    write(fid,'(A)') 'ecpbasis '//TRIM(elem(i))//' '//TRIM(elem(i))//'.gen'
   else
    write(fid,'(A)') 'basis  '//TRIM(elem(i))//' '//TRIM(elem(i))//'.gen'
   end if
  else ! not allocated
    write(fid,'(A)') 'basis  '//TRIM(elem(i))//' '//TRIM(elem(i))//'.gen'
  end if
 end do ! for i

 write(fid,'(A)') 'print  normal'
 write(fid,'(A)') 'end'

 write(fid,'(/,A)') '$doqm'

 write(fid,'(/,A)') '$xml'
 write(fid,'(A)') 'ow '//TRIM(xmlname)
 write(fid,'(A)') 'denmatw'
 write(fid,'(A)') 'cw'
 write(fid,'(A)') 'end'
 close(fid)
 return
end subroutine write_qm4d_inp

! write/create basis files (and ECP, if any) for each element
subroutine write_basis_files_for_each_elem()
 use fch_content, only: natom, elem, ncontr, prim_exp, contr_coeff, contr_coeff_sp,&
  shell_type, shell2atom_map, prim_per_shell, LPSkip, LMax, RNFroz, KFirst, &
  KLast, NLP, CLP, ZLP
 implicit none
 integer :: i, j, k0, k, m, n, p, q, fid
 character(len=1), parameter :: am_type(0:6)=['S','P','D','F','G','H','I']
 character(len=2) :: str
 character(len=30) :: basname

 ! print basis sets into the .inp file
 k = 0 ! current index in array shell_type
 p = 0 ! current index in array prim_exp

 do i = 1, natom, 1
  if(i > 1) then
   if(ANY(elem(1:i-1) == elem(i))) then
    if(i == natom) exit
    k0 = k
    do j = k0+1, ncontr, 1
     if(shell2atom_map(j) == i+1) exit
     k = k + 1
     p = p + prim_per_shell(j)
    end do ! for j
    cycle
   end if
  end if

  basname = TRIM(elem(i))//'.gen'
  open(newunit=fid,file=TRIM(basname),status='replace')
  write(fid,'(A)') '****'
  write(fid,'(A)') elem(i)//' 0'

  k = k + 1
  k0 = k
  do j = k0, ncontr, 1
   if(shell2atom_map(j) == i+1) exit
   m = shell_type(j); n = prim_per_shell(j)
   if(m == -1) then
    str = 'SP'
   else
    str = am_type(m)//' '
   end if
   write(fid,'(A2,1X,I2,A)') str, n, ' 1.00'

   do q = p+1, p+n, 1
    if(allocated(contr_coeff_sp)) then
     if(DABS(contr_coeff_sp(q)) > 1e-6) then
      write(fid,'(3(2X,ES15.8))') prim_exp(q), contr_coeff(q), contr_coeff_sp(q)
     else
      write(fid,'(2(2X,ES15.8))') prim_exp(q), contr_coeff(q)
     end if
    else ! no 'SP'
     write(fid,'(2(2X,ES15.8))') prim_exp(q), contr_coeff(q)
    end if
   end do ! for q

   p = p + n ! update p
  end do ! for j

  write(fid,'(A)') '****'
  k = j - 1 ! update k

  if(allocated(LPSkip)) then
   if(LPSkip(i) == 0) call prt_ecp_in_gau_format(fid, i)
  end if
  close(fid)
 end do ! for i

 deallocate(elem, shell2atom_map, prim_per_shell, prim_exp, contr_coeff)
 if(allocated(contr_coeff_sp)) deallocate(contr_coeff_sp)
 if(allocated(LPSkip)) deallocate(LPSkip)
 if(allocated(LMax)) deallocate(LMax)
 if(allocated(RNFroz)) deallocate(RNFroz)
 if(allocated(KFirst)) deallocate(KFirst)
 if(allocated(KLast)) deallocate(KLast)
 if(allocated(NLP)) deallocate(NLP)
 if(allocated(CLP)) deallocate(CLP)
 if(allocated(ZLP)) deallocate(ZLP)
 return
end subroutine write_basis_files_for_each_elem

! print ECP data in Gaussian format
subroutine prt_ecp_in_gau_format(fid, i)
 use fch_content, only: elem, LMax, RNFroz, KFirst, KLast, NLP, CLP, ZLP
 implicit none
 integer :: j, n, n1, n2
 integer, intent(in) :: fid, i
 character(len=1) :: str
 character(len=1), parameter :: am_type1(0:6)=['s','p','d','f','g','h','i']

 write(fid,'(/,A)') elem(i)//' 0'
 write(fid,'(A,I3,1X,I4)') TRIM(elem(i))//'-ECP ',LMax(i),INT(RNFroz(i))
 str = am_type1(LMax(i))

 do j = 1, 10, 1
  n1 = KFirst(i,j); n2 = KLast(i,j)
  if(n1 == 0) exit
  if(j == 1) then
   write(fid,'(A)') str//' potential'
  else
   write(fid,'(A)') am_type1(j-2)//'-'//str//' potential'
  end if
  write(fid,'(I3)') n2-n1+1
  do n = n1, n2, 1
   write(fid,'(I0,2(3X,ES15.8))') NLP(n), ZLP(n), CLP(n)
  end do ! for n
 end do ! for j

 return
end subroutine prt_ecp_in_gau_format

subroutine fch2qm4d_permute_6d(nif,coeff)
 implicit none
 integer :: i
 integer, parameter :: order(6) = [1, 4, 5, 2, 6, 3]
 integer, intent(in) :: nif
 real(kind=8), intent(inout) :: coeff(6,nif)
 real(kind=8), allocatable :: coeff2(:,:)
! From: the order of Cartesian d functions in Gaussian
! To: the order of Cartesian d functions in QM4D
! 1  2  3  4  5  6
! XX,YY,ZZ,XY,XZ,YZ
! XX,XY,XZ,YY,YZ,ZZ

 allocate(coeff2(6,nif), source=0d0)
 forall(i = 1:6) coeff2(i,:) = coeff(order(i),:)
 coeff = coeff2
 deallocate(coeff2)
 return
end subroutine fch2qm4d_permute_6d

subroutine fch2qm4d_permute_10f(nif,coeff)
 implicit none
 integer :: i
 integer, parameter :: order(10) = [1, 5, 6, 4, 10, 7, 2, 9, 8, 3]
 integer, intent(in) :: nif
 real(kind=8), intent(inout) :: coeff(10,nif)
 real(kind=8), allocatable :: coeff2(:,:)
! From: the order of Cartesian f functions in Gaussian
! To: the order of Cartesian f functions in QM4D
! 1   2   3   4   5   6   7   8   9   10
! XXX,YYY,ZZZ,XYY,XXY,XXZ,XZZ,YZZ,YYZ,XYZ
! XXX,XXY,XXZ,XYY,XYZ,XZZ,YYY,YYZ,YZZ,ZZZ

 allocate(coeff2(10,nif), source=0d0)
 forall(i = 1:10) coeff2(i,:) = coeff(order(i),:)
 coeff = coeff2
 deallocate(coeff2)
 return
end subroutine fch2qm4d_permute_10f


text:
      subroutine VV_HZZ_v(p,msqv)
      implicit none
************************************************************************
*     Author: J. M. Campbell                                           *
*     June, 2002.                                                      *
*                                                                      *
*     Weak Boson Fusion : sums up WW and ZZ contributions              *
*     This routine calculates the virtual matrix element squared       *
*     for the process:                                                 *
*     q(-p1) + q(-p2) --> H(p34) + q(p5) + q(p6)                       *
*                           |                                          *
*                           |                                          *
*                           |                                          *
*                           ---> W+(nu(p3)+e+(p4))+W-(e-(p5)+nub(p6))  *
*                                                                      *
************************************************************************
      include 'constants.f'
      include 'qcdcouple.f'
      include 'epinv.f'
      include 'epinv2.f'
      include 'scale.f'
      include 'scheme.f'
      double precision dot,p(mxpart,4),
     . msq0(-nf:nf,-nf:nf),msqv(-nf:nf,-nf:nf),xl17,xl28,facv
      integer j,k

c--- this result is in the MS-bar scheme
      scheme='dred'

* Virtual matrix elements are simply lowest order multiplied by factor
* Need to check this for pisq (cf. qqb_w_v.f) and I think there should
* be an additional +ason2pi*cf*(1+1) from DRED compared to MSBAR
      call VV_HZZ(p,msq0)
      xl17=dlog(-2d0*dot(p,1,7)/musq)
      xl28=dlog(-2d0*dot(p,2,8)/musq)
      facv=ason2pi*cf*(
     . -4d0*EPINV*EPINV2-(6d0-2d0*(xl17+xl28))*EPINV
     . +3d0*(xl17+xl28)-(xl17**2+xl28**2)-14d0)

      do j=-nf,nf
      do k=-nf,nf
        msqv(j,k)=facv*msq0(j,k)
      enddo
      enddo

      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
C
C	$Id: csdelnod.f,v 1.5 2008-07-27 03:10:06 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE CSDELNOD (K, N,X,Y,Z,LIST,LPTR,LEND,LNEW,LWK,
     .                   IWK, IER)
      INTEGER K, N, LIST(*), LPTR(*), LEND(*), LNEW, LWK,
     .        IWK(2,*), IER
      DOUBLE PRECISION X(*), Y(*), Z(*)
C
C***********************************************************
C
C                                              From STRIPACK
C                                            Robert J. Renka
C                                  Dept. of Computer Science
C                                       Univ. of North Texas
C                                           renka@cs.unt.edu
C                                                   11/30/99
C
C   This subroutine deletes node K (along with all arcs
C incident on node K) from a triangulation of N nodes on the
C unit sphere, and inserts arcs as necessary to produce a
C triangulation of the remaining N-1 nodes.  If a Delaunay
C triangulation is input, a Delaunay triangulation will
C result, and thus, CSDELNOD reverses the effect of a call to
C Subroutine CSADDNOD.
C
C
C On input:
C
C       K = Index (for X, Y, and Z) of the node to be
C           deleted.  1 .LE. K .LE. N.
C
C K is not altered by this routine.
C
C       N = Number of nodes in the triangulation on input.
C           N .GE. 4.  Note that N will be decremented
C           following the deletion.
C
C       X,Y,Z = Arrays of length N containing the Cartesian
C               coordinates of the nodes in the triangula-
C               tion.
C
C       LIST,LPTR,LEND,LNEW = Data structure defining the
C                             triangulation.  Refer to Sub-
C                             routine CSTRMESH.
C
C       LWK = Number of columns reserved for IWK.  LWK must
C             be at least NNB-3, where NNB is the number of
C             neighbors of node K, including an extra
C             pseudo-node if K is a boundary node.
C
C       IWK = Integer work array dimensioned 2 by LWK (or
C             array of length .GE. 2*LWK).
C
C On output:
C
C       N = Number of nodes in the triangulation on output.
C           The input value is decremented unless 1 .LE. IER
C           .LE. 4.
C
C       X,Y,Z = Updated arrays containing nodal coordinates
C               (with elements K+1,...,N+1 shifted up one
C               position, thus overwriting element K) unless
C               1 .LE. IER .LE. 4.
C
C       LIST,LPTR,LEND,LNEW = Updated triangulation data
C                             structure reflecting the dele-
C                             tion unless 1 .LE. IER .LE. 4.
C                             Note that the data structure
C                             may have been altered if IER >
C                             3.
C
C       LWK = Number of IWK columns required unless IER = 1
C             or IER = 3.
C
C       IWK = Indexes of the endpoints of the new arcs added
C             unless LWK = 0 or 1 .LE. IER .LE. 4.  (Arcs
C             are associated with columns, or pairs of
C             adjacent elements if IWK is declared as a
C             singly-subscripted array.)
C
C       IER = Error indicator:
C             IER = 0 if no errors were encountered.
C             IER = 1 if K or N is outside its valid range
C                     or LWK < 0 on input.
C             IER = 2 if more space is required in IWK.
C                     Refer to LWK.
C             IER = 3 if the triangulation data structure is
C                     invalid on input.
C             IER = 4 if K indexes an interior node with
C                     four or more neighbors, none of which
C                     can be swapped out due to collineari-
C                     ty, and K cannot therefore be deleted.
C             IER = 5 if an error flag (other than IER = 1)
C                     was returned by CSOPTIM.  An error
C                     message is written to the standard
C                     output unit in this case.
C             IER = 6 if error flag 1 was returned by CSOPTIM.
C                     This is not necessarily an error, but
C                     the arcs may not be optimal.
C
C   Note that the deletion may result in all remaining nodes
C being collinear.  This situation is not flagged.
C
C Modules required by CSDELNOD:  CSDELNB, CSLEFT, CSLSTPTR, CSNBCNT,
C                                CSOPTIM, CSSWAP, CSSWPTST
C
C Intrinsic function called by CSDELNOD:  ABS
C
C***********************************************************
C
      INTEGER CSLSTPTR, CSNBCNT
      INTEGER I, IERR, IWL, J, LNW, LP, LP21, LPF, LPH, LPL,
     .        LPL2, LPN, LWKL, N1, N2, NFRST, NIT, NL, NN,
     .        NNB, NR
      LOGICAL CSLEFT
      LOGICAL BDRY
      DOUBLE PRECISION X1, X2, XL, XR, Y1, Y2, YL, YR, Z1,
     .                 Z2, ZL, ZR
C
C Local parameters:
C
C BDRY =    Logical variable with value TRUE iff N1 is a
C             boundary node
C I,J =     DO-loop indexes
C IERR =    Error flag returned by CSOPTIM
C IWL =     Number of IWK columns containing arcs
C LNW =     Local copy of LNEW
C LP =      LIST pointer
C LP21 =    LIST pointer returned by CSSWAP
C LPF,LPL = Pointers to the first and last neighbors of N1
C LPH =     Pointer (or flag) returned by CSDELNB
C LPL2 =    Pointer to the last neighbor of N2
C LPN =     Pointer to a neighbor of N1
C LWKL =    Input value of LWK
C N1 =      Local copy of K
C N2 =      Neighbor of N1
C NFRST =   First neighbor of N1:  LIST(LPF)
C NIT =     Number of iterations in CSOPTIM
C NR,NL =   Neighbors of N1 preceding (to the right of) and
C             following (to the left of) N2, respectively
C NN =      Number of nodes in the triangulation
C NNB =     Number of neighbors of N1 (including a pseudo-
C             node representing the boundary if N1 is a
C             boundary node)
C X1,Y1,Z1 = Coordinates of N1
C X2,Y2,Z2 = Coordinates of N2
C XL,YL,ZL = Coordinates of NL
C XR,YR,ZR = Coordinates of NR
C
C
C Set N1 to K and NNB to the number of neighbors of N1 (plus
C   one if N1 is a boundary node), and test for errors.  LPF
C   and LPL are LIST indexes of the first and last neighbors
C   of N1, IWL is the number of IWK columns containing arcs,
C   and BDRY is TRUE iff N1 is a boundary node.
C
      N1 = K
      NN = N
      IF (N1 .LT. 1  .OR.  N1 .GT. NN  .OR.  NN .LT. 4  .OR.
     .    LWK .LT. 0) GO TO 21
      LPL = LEND(N1)
      LPF = LPTR(LPL)
      NNB = CSNBCNT(LPL,LPTR)
      BDRY = LIST(LPL) .LT. 0
      IF (BDRY) NNB = NNB + 1
      IF (NNB .LT. 3) GO TO 23
      LWKL = LWK
      LWK = NNB - 3
      IF (LWKL .LT. LWK) GO TO 22
      IWL = 0
      IF (NNB .EQ. 3) GO TO 3
C
C Initialize for loop on arcs N1-N2 for neighbors N2 of N1,
C   beginning with the second neighbor.  NR and NL are the
C   neighbors preceding and following N2, respectively, and
C   LP indexes NL.  The loop is exited when all possible
C   swaps have been applied to arcs incident on N1.
C
      X1 = X(N1)
      Y1 = Y(N1)
      Z1 = Z(N1)
      NFRST = LIST(LPF)
      NR = NFRST
      XR = X(NR)
      YR = Y(NR)
      ZR = Z(NR)
      LP = LPTR(LPF)
      N2 = LIST(LP)
      X2 = X(N2)
      Y2 = Y(N2)
      Z2 = Z(N2)
      LP = LPTR(LP)
C
C Top of loop:  set NL to the neighbor following N2.
C
    1 NL = ABS(LIST(LP))
      IF (NL .EQ. NFRST  .AND.  BDRY) GO TO 3
      XL = X(NL)
      YL = Y(NL)
      ZL = Z(NL)
C
C   Test for a convex quadrilateral.  To avoid an incorrect
C     test caused by collinearity, use the fact that if N1
C     is a boundary node, then N1 CSLEFT NR->NL and if N2 is
C     a boundary node, then N2 CSLEFT NL->NR.
C
      LPL2 = LEND(N2)
      IF ( .NOT. ((BDRY  .OR.  CSLEFT(XR,YR,ZR,XL,YL,ZL,X1,Y1,
     .      Z1))  .AND.  (LIST(LPL2) .LT. 0  .OR.
     .      CSLEFT(XL,YL,ZL,XR,YR,ZR,X2,Y2,Z2))) ) THEN
C
C   Nonconvex quadrilateral -- no swap is possible.
C
        NR = N2
        XR = X2
        YR = Y2
        ZR = Z2
        GO TO 2
      ENDIF
C
C   The quadrilateral defined by adjacent triangles
C     (N1,N2,NL) and (N2,N1,NR) is convex.  Swap in
C     NL-NR and store it in IWK unless NL and NR are
C     already adjacent, in which case the swap is not
C     possible.  Indexes larger than N1 must be decremented
C     since N1 will be deleted from X, Y, and Z.
C
      CALL CSSWAP (NL,NR,N1,N2, LIST,LPTR,LEND, LP21)
      IF (LP21 .EQ. 0) THEN
        NR = N2
        XR = X2
        YR = Y2
        ZR = Z2
        GO TO 2
      ENDIF
      IWL = IWL + 1
      IF (NL .LE. N1) THEN
        IWK(1,IWL) = NL
      ELSE
        IWK(1,IWL) = NL - 1
      ENDIF
      IF (NR .LE. N1) THEN
        IWK(2,IWL) = NR
      ELSE
        IWK(2,IWL) = NR - 1
      ENDIF
C
C   Recompute the LIST indexes and NFRST, and decrement NNB.
C
      LPL = LEND(N1)
      NNB = NNB - 1
      IF (NNB .EQ. 3) GO TO 3
      LPF = LPTR(LPL)
      NFRST = LIST(LPF)
      LP = CSLSTPTR(LPL,NL,LIST,LPTR)
      IF (NR .EQ. NFRST) GO TO 2
C
C   NR is not the first neighbor of N1.
C     Back up and test N1-NR for a swap again:  Set N2 to
C     NR and NR to the previous neighbor of N1 -- the
C     neighbor of NR which follows N1.  LP21 points to NL
C     as a neighbor of NR.
C
      N2 = NR
      X2 = XR
      Y2 = YR
      Z2 = ZR
      LP21 = LPTR(LP21)
      LP21 = LPTR(LP21)
      NR = ABS(LIST(LP21))
      XR = X(NR)
      YR = Y(NR)
      ZR = Z(NR)
      GO TO 1
C
C   Bottom of loop -- test for termination of loop.
C
    2 IF (N2 .EQ. NFRST) GO TO 3
      N2 = NL
      X2 = XL
      Y2 = YL
      Z2 = ZL
      LP = LPTR(LP)
      GO TO 1
C
C Delete N1 and all its incident arcs.  If N1 is an interior
C   node and either NNB > 3 or NNB = 3 and N2 CSLEFT NR->NL,
C   then N1 must be separated from its neighbors by a plane
C   containing the origin -- its removal reverses the effect
C   of a call to CSCOVSPH, and all its neighbors become
C   boundary nodes.  This is achieved by treating it as if
C   it were a boundary node (setting BDRY to TRUE, changing
C   a sign in LIST, and incrementing NNB).
C
    3 IF (.NOT. BDRY) THEN
        IF (NNB .GT. 3) THEN
          BDRY = .TRUE.
        ELSE
          LPF = LPTR(LPL)
          NR = LIST(LPF)
          LP = LPTR(LPF)
          N2 = LIST(LP)
          NL = LIST(LPL)
          BDRY = CSLEFT(X(NR),Y(NR),Z(NR),X(NL),Y(NL),Z(NL),
     .                X(N2),Y(N2),Z(N2))
        ENDIF
        IF (BDRY) THEN
C
C   IF a boundary node already exists, then N1 and its
C     neighbors cannot be converted to boundary nodes.
C     (They must be collinear.)  This is a problem if
C     NNB > 3.
C
          DO 4 I = 1,NN
            IF (LIST(LEND(I)) .LT. 0) THEN
              BDRY = .FALSE.
              GO TO 5
            ENDIF
    4       CONTINUE
          LIST(LPL) = -LIST(LPL)
          NNB = NNB + 1
        ENDIF
      ENDIF
    5 IF (.NOT. BDRY  .AND.  NNB .GT. 3) GO TO 24
C
C Initialize for loop on neighbors.  LPL points to the last
C   neighbor of N1.  LNEW is stored in local variable LNW.
C
      LP = LPL
      LNW = LNEW
C
C Loop on neighbors N2 of N1, beginning with the first.
C
    6 LP = LPTR(LP)
        N2 = ABS(LIST(LP))
        CALL CSDELNB (N2,N1,N, LIST,LPTR,LEND,LNW, LPH)
        IF (LPH .LT. 0) GO TO 23
C
C   LP and LPL may require alteration.
C
        IF (LPL .EQ. LNW) LPL = LPH
        IF (LP .EQ. LNW) LP = LPH
        IF (LP .NE. LPL) GO TO 6
C
C Delete N1 from X, Y, Z, and LEND, and remove its adjacency
C   list from LIST and LPTR.  LIST entries (nodal indexes)
C   which are larger than N1 must be decremented.
C
      NN = NN - 1
      IF (N1 .GT. NN) GO TO 9
      DO 7 I = N1,NN
        X(I) = X(I+1)
        Y(I) = Y(I+1)
        Z(I) = Z(I+1)
        LEND(I) = LEND(I+1)
    7   CONTINUE
C
      DO 8 I = 1,LNW-1
        IF (LIST(I) .GT. N1) LIST(I) = LIST(I) - 1
        IF (LIST(I) .LT. -N1) LIST(I) = LIST(I) + 1
    8   CONTINUE
C
C   For LPN = first to last neighbors of N1, delete the
C     preceding neighbor (indexed by LP).
C
C   Each empty LIST,LPTR location LP is filled in with the
C     values at LNW-1, and LNW is decremented.  All pointers
C     (including those in LPTR and LEND) with value LNW-1
C     must be changed to LP.
C
C  LPL points to the last neighbor of N1.
C
    9 IF (BDRY) NNB = NNB - 1
      LPN = LPL
      DO 13 J = 1,NNB
        LNW = LNW - 1
        LP = LPN
        LPN = LPTR(LP)
        LIST(LP) = LIST(LNW)
        LPTR(LP) = LPTR(LNW)
        IF (LPTR(LPN) .EQ. LNW) LPTR(LPN) = LP
        IF (LPN .EQ. LNW) LPN = LP
        DO 10 I = NN,1,-1
          IF (LEND(I) .EQ. LNW) THEN
            LEND(I) = LP
            GO TO 11
          ENDIF
   10     CONTINUE
C
   11   DO 12 I = LNW-1,1,-1
          IF (LPTR(I) .EQ. LNW) LPTR(I) = LP
   12     CONTINUE
   13   CONTINUE
C
C Update N and LNEW, and optimize the patch of triangles
C   containing K (on input) by applying swaps to the arcs
C   in IWK.
C
      N = NN
      LNEW = LNW
      IF (IWL .GT. 0) THEN
        NIT = 4*IWL
        CALL CSOPTIM (X,Y,Z,IWL, LIST,LPTR,LEND,NIT,IWK, IERR)
        IF (IERR .NE. 0  .AND.  IERR .NE. 1) GO TO 25
        IF (IERR .EQ. 1) GO TO 26
      ENDIF
C
C Successful termination.
C
      IER = 0
      RETURN
C
C Invalid input parameter.
C
   21 IER = 1
      RETURN
C
C Insufficient space reserved for IWK.
C
   22 IER = 2
      RETURN
C
C Invalid triangulation data structure.  NNB < 3 on input or
C   N2 is a neighbor of N1 but N1 is not a neighbor of N2.
C
   23 IER = 3
      RETURN
C
C N1 is interior but NNB could not be reduced to 3.
C
   24 IER = 4
      RETURN
C
C Error flag (other than 1) returned by CSOPTIM.
C
   25 IER = 5
      WRITE (*,100) NIT, IERR
  100 FORMAT (//5X,'*** Error in CSOPTIM (called from ',
     .        'CSDELNOD):  NIT = ',I4,', IER = ',I1,' ***'/)
      RETURN
C
C Error flag 1 returned by CSOPTIM.
C
   26 IER = 6
      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
!
! NAME:
!       MWLBL_Liebe93
!
! PURPOSE:
!       Module containing routines to calculate microwave atmospheric
!       attenuation according to the Liebe 93 MPM model.
!
!
! CREATION HISTORY:
!       This module has been translated from FORTRAN-77 source code
!       and include files written by:
!         L. Phalippou ECMWF 13-Dec-1993
!         Peter Rayer UKMO 19-Dec-2000
!         Roger Saunders UKMO
!
!       Written by:     Paul van Delst, 21-Apr-2002
!                       paul.vandelst@noaa.gov
!

MODULE MWLBL_Liebe93

  ! -----------------
  ! Environment setup
  ! -----------------
  ! Module usage
  USE Type_Kinds     , ONLY: fp
  USE Message_Handler, ONLY: SUCCESS, FAILURE, WARNING, Display_Message
  USE Liebe92_Coefficients
!  USE Liebe93_Coefficients
  ! Disable all implicit typing
  IMPLICIT NONE

  ! ------------
  ! Visibilities
  ! ------------
  PRIVATE
  PUBLIC :: Liebe93

  ! -------------------
  ! Procedure overloads
  ! -------------------
  INTERFACE Liebe93
    MODULE PROCEDURE Liebe93_By_Layer
    MODULE PROCEDURE Liebe93_By_Frequency
  END INTERFACE Liebe93

  ! -----------------
  ! Module parameters
  ! -----------------
  CHARACTER(*), PARAMETER :: MODULE_VERSION_ID = &
  ! Literal constants
  REAL(fp), PARAMETER :: ZERO           = 0.0_fp
  REAL(fp), PARAMETER :: ZEROpointEIGHT = 0.8_fp
  REAL(fp), PARAMETER :: ONE            = 1.0_fp
  REAL(fp), PARAMETER :: ONEpointONE    = 1.1_fp
  REAL(fp), PARAMETER :: ONEpointFIVE   = 1.5_fp
  REAL(fp), PARAMETER :: THREE          = 3.0_fp
  REAL(fp), PARAMETER :: THREEpointFIVE = 3.5_fp
  REAL(fp), PARAMETER :: FOURpointFIVE  = 4.5_fp
  REAL(fp), PARAMETER :: SEVENpointFIVE = 7.5_fp
  REAL(fp), PARAMETER :: TOLERANCE      = EPSILON(ZERO)
  ! Conversion factor for attenuation computation.
  REAL(fp), PARAMETER :: ATT_K = 0.1820_fp
  ! Min/max temperatures in Kelvin
  REAL(fp), PARAMETER :: MIN_TEMPERATURE = 150.0_fp
  REAL(fp), PARAMETER :: MAX_TEMPERATURE = 350.0_fp
  ! Min/max frequencies in GHz
  REAL(fp), PARAMETER :: MIN_FREQUENCY = ONE
  REAL(fp), PARAMETER :: MAX_FREQUENCY = 1000.0_fp
  ! Default message string length
  INTEGER,  PARAMETER :: ML = 256

CONTAINS

!##################################################################################
!##################################################################################
!##                                                                              ##
!##                           ## PUBLIC MODULE ROUTINES ##                       ##
!##                                                                              ##
!##################################################################################
!##################################################################################

!------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       Liebe93
!
! PURPOSE:
!       Function to calculate the atmospheric attentuation in the microwave
!       according to Liebe's MPM93 model.
!
! CALLING SEQUENCE:
!       Error_Status = Liebe93( &
!         Frequency       , &  ! Input
!         Dry_Air_Pressure, &  ! Input
!         H2O_Pressure    , &  ! Input
!         Temperature     , &  ! Input
!         Quiet                    = Quiet                   , &  ! Optional input
!         WetLine_Attenuation      = WetLine_Attenuation     , &  ! Optional output
!         WetContinuum_Attenuation = WetContinuum_Attenuation, &  ! Optional output
!         DryLine_Attenuation      = DryLine_Attenuation     , &  ! Optional output
!         DryContinuum_Attenuation = DryContinuum_Attenuation  )  ! Optional output
!
!
!       Two forms of argument specification are available - by FREQUENCY
!       or by LAYER.
!
!       If by FREQUENCY, then multiple frequencies are processed for
!       a single layer:
!         Frequency        - Rank-1, size n_Frequencies
!         Dry_Air_Pressure - Scalar
!         H2O_Pressure     - Scalar
!         Temperature      - Scalar
!         Attenuation      - Rank-1, size n_Frequencies
!
!       If by LAYER, then a single frequency is processed for multiple
!       layers:
!         Frequency        - Scalar
!         Dry_Air_Pressure - Rank-1, size n_Layers
!         H2O_Pressure     - Rank-1, size n_Layers
!         Temperature      - Rank-1, size n_Layers
!         Attenuation      - Rank-1, size n_Layers
!
!
! INPUTS:
!       Frequency:                  Frequency for which the atmospheric attenuation
!                                   is required
!                                   UNITS:      GHz
!                                   TYPE:       REAL(fp)
!                                   DIMENSION:  Rank-1  (by FREQUENCY)
!                                                 OR
!                                               Scalar  (by LAYER)
!                                   ATTRIBUTES: INTENT(IN)
!
!       Dry_Air_Pressure:           Dry air partial pressure profile.
!                                   UNITS:      hectoPascals, hPa
!                                   TYPE:       REAL(fp)
!                                   DIMENSION:  Scalar  (by FREQUENCY)
!                                                 OR
!                                               Rank-1  (by LAYER)
!                                   ATTRIBUTES: INTENT(IN)
!
!       H2O_Pressure:               Water vapor partial pressure profile.
!                                   UNITS:      hectoPascals, hPa
!                                   TYPE:       REAL(fp)
!                                   DIMENSION:  Same as DRY_AIR_PRESSURE argument
!                                   ATTRIBUTES: INTENT(IN)
!
!       Temperature:                Temperature profile.
!                                   UNITS:      Kelvin
!                                   TYPE:       REAL(fp)
!                                   DIMENSION:  Same as DRY_AIR_PRESSURE argument
!                                   ATTRIBUTES: INTENT(IN)
!
! OPTIONAL INPUTS:
!       Quiet:                      Set this logical argument to suppress WARNING and
!                                   INFORMATION messages being printed to stdout
!                                   If == .FALSE., messages are OUTPUT [DEFAULT].
!                                      == .TRUE.,  messages are SUPPRESSED.
!                                   If not specified, default is .FALSE.
!                                   UNITS:      N/A
!                                   TYPE:       LOGICAL
!                                   DIMENSION:  Scalar
!                                   ATTRIBUTES: INTENT(IN), OPTIONAL
!
! OPTIONAL OUTPUTS:
!       WetLine_Attenuation:        Power attenuation due to absorption by water
!                                   vapour absorption lines.
!                                   UNITS:      dB/km
!                                   TYPE:       REAL(fp)
!                                   DIMENSION:  Rank-1  (n_Frequencies, by FREQUENCY)
!                                                 OR
!                                               Rank-1  (n_Layers, by LAYER)
!                                   ATTRIBUTES: INTENT(OUT), OPTIONAL
!
!       WetContinuum_Attenuation:   Power attenuation due to water vapour continuum
!                                   absorption.
!                                   UNITS:      dB/km
!                                   TYPE:       REAL(fp)
!                                   DIMENSION:  Same as WETLINE_ATTENUATION argument.
!                                   ATTRIBUTES: INTENT(OUT), OPTIONAL
!
!       DryLine_Attenuation:        Power attenuation due to absorption by oxygen
!                                   absorption lines.
!                                   UNITS:      dB/km
!                                   TYPE:       REAL(fp)
!                                   DIMENSION:  Same as WETLINE_ATTENUATION argument.
!                                   ATTRIBUTES: INTENT(OUT), OPTIONAL
!
!       DryContinuum_Attenuation:   Power attenuation due to non-resonant dry air
!                                   absorption.
!                                   UNITS:      dB/km
!                                   TYPE:       REAL(fp)
!                                   DIMENSION:  Same as WETLINE_ATTENUATION argument.
!                                   ATTRIBUTES: INTENT(OUT), OPTIONAL
!
! FUNCTION RESULT:
!       Error_Status:               The return value is an integer defining the error status.
!                                   status. The error codes are defined in the Message_Handler
!                                   module.
!                                   If == SUCCESS the calculation was successful
!                                      == FAILURE an unrecoverable error occurred.
!                                   UNITS:      N/A
!                                   TYPE:       INTEGER
!                                   DIMENSION:  Scalar
!
! PROCEDURE:
!       Liebe, H.J., G.A. Hufford, and M.G. Cotton, 1993: "Propagation
!         modelling of moist air and suspended water/ice particles at
!         frequencies below 1000 GHz". In AGARD Conference Proceedings
!         542, Atmospheric propagation effects through natural and man-made
!         obscurants for visible through MM-wave radiation, pp 3.1-3.10.
!         Presented at the Electromagnetic Wave Propagation Panel Symposium,
!         Palma de Mallorca, Spain, 17-20 May 1993.
!
! CREATION HISTORY:
!       F77 version written by : L. Phalippou ECMWF 13-Dec-1993
!       F77 version modified by: Peter Rayer UKMO 19-Dec-2000
!                                Roger Saunders UKMO
!
!       F90 version translation: Paul van Delst, CIMSS/SSEC 21-Apr-2002
!                                paul.vandelst@ssec.wisc.edu
!
!------------------------------------------------------------------------------

  ! --------------------------------------
  ! Multi-Frequency, Single-Layer function
  ! --------------------------------------
  FUNCTION Liebe93_By_Frequency( &
    Frequency               , &  ! Input
    Dry_Air_Pressure        , &  ! Input
    H2O_Pressure            , &  ! Input
    Temperature             , &  ! Input
    Quiet,                    &  ! Optional input
    WetLine_Attenuation     , &  ! Optional output
    WetContinuum_Attenuation, &  ! Optional output
    DryLine_Attenuation     , &  ! Optional output
    DryContinuum_Attenuation) &  ! Optional output
  RESULT( err_stat )
    ! Arguments
    REAL(fp),           INTENT(IN)  :: Frequency(:)
    REAL(fp),           INTENT(IN)  :: Dry_Air_Pressure
    REAL(fp),           INTENT(IN)  :: H2O_Pressure
    REAL(fp),           INTENT(IN)  :: Temperature
    LOGICAL , OPTIONAL, INTENT(IN)  :: Quiet
    REAL(fp), OPTIONAL, INTENT(OUT) :: WetLine_Attenuation(:)
    REAL(fp), OPTIONAL, INTENT(OUT) :: WetContinuum_Attenuation(:)
    REAL(fp), OPTIONAL, INTENT(OUT) :: DryLine_Attenuation(:)
    REAL(fp), OPTIONAL, INTENT(OUT) :: DryContinuum_Attenuation(:)
    ! Function result
    INTEGER :: err_stat
    ! Local parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'Liebe93::By_Frequency'
    ! Local variables
    CHARACTER(ML) :: msg
    LOGICAL :: noisy
    LOGICAL :: Compute_WetLine
    LOGICAL :: Compute_WetContinuum
    LOGICAL :: Compute_DryLine
    LOGICAL :: Compute_DryContinuum
    INTEGER :: n_Frequencies, l
    REAL(fp) :: Theta
    REAL(fp) :: O2_s1, Common_s1, O2_g1, O2_d1, H2O_s1
    REAL(fp) :: H2O_Line_Absorption
    REAL(fp) :: WetContinuum_Absorption
    REAL(fp) :: O2_Line_Absorption
    REAL(fp) :: DryContinuum_Absorption

    ! Set up
    err_stat = SUCCESS
    ! ...Process keywords
    noisy = .TRUE.
    IF ( PRESENT(Quiet) ) noisy = .NOT.Quiet

    ! Check output argument validity
    n_Frequencies = SIZE(Frequency)
    ! ...The water vapour line attenuation
    Compute_WetLine = .FALSE.
    IF ( PRESENT(WetLine_Attenuation) ) THEN
      Compute_WetLine = .TRUE.
      IF ( SIZE(WetLine_Attenuation) /= n_Frequencies ) THEN
        err_stat = FAILURE
        msg = 'WetLine_Attenuation argument array has inconsistent size.'
        CALL Display_Message( ROUTINE_NAME, msg, err_stat ); RETURN
      END IF
    END IF
    ! ...The water vapour continuum attenuation
    Compute_WetContinuum = .FALSE.
    IF ( PRESENT(WetContinuum_Attenuation) ) THEN
      Compute_WetContinuum = .TRUE.
      IF ( SIZE(WetContinuum_Attenuation) /= n_Frequencies ) THEN
        err_stat = FAILURE
        msg = 'WetContinuum_Attenuation argument array has inconsistent size.'
        CALL Display_Message( ROUTINE_NAME, msg, err_stat ); RETURN
      END IF
    END IF
    ! ...The water vapour line attenuation
    Compute_DryLine = .FALSE.
    IF ( PRESENT(DryLine_Attenuation) ) THEN
      Compute_DryLine = .TRUE.
      IF ( SIZE(DryLine_Attenuation) /= n_Frequencies ) THEN
        err_stat = FAILURE
        msg = 'DryLine_Attenuation argument array has inconsistent size.'
        CALL Display_Message( ROUTINE_NAME, msg, err_stat ); RETURN
      END IF
    END IF
    ! ...The water vapour continuum attenuation
    Compute_DryContinuum = .FALSE.
    IF ( PRESENT( DryContinuum_Attenuation ) ) THEN
      Compute_DryContinuum = .TRUE.
      IF ( SIZE( DryContinuum_Attenuation ) /= n_Frequencies ) THEN
        err_stat = FAILURE
        msg = 'DryContinuum_Attenuation argument array has inconsistent size.'
        CALL Display_Message( ROUTINE_NAME, msg, err_stat ); RETURN
      END IF
    END IF
    ! ...If no output arguments at all, return
    IF ( ( .NOT. Compute_WetLine      ) .AND. &
         ( .NOT. Compute_WetContinuum ) .AND. &
         ( .NOT. Compute_DryLine      ) .AND. &
         ( .NOT. Compute_DryContinuum )       ) RETURN
    ! ...Check for invalid frequencies
    IF ( ANY( Frequency < ZERO ) ) THEN
      err_stat = FAILURE
        msg = 'Input frequencies must be > 0.0.'
        CALL Display_Message( ROUTINE_NAME, msg, err_stat ); RETURN
    END IF
    ! ...Check for invalid profile data
    IF ( Dry_Air_Pressure < TOLERANCE .OR. &
         H2O_Pressure     < TOLERANCE .OR. &
         Temperature      < TOLERANCE      ) THEN
      err_stat = FAILURE
      msg = 'Input tressure/temperature arguments must be > 0.0.'
      CALL Display_Message( ROUTINE_NAME, msg, err_stat ); RETURN
    END IF

    ! Issue warnings for data outside range
    IF ( noisy )  THEN

      ! Frequencies < min value or > max_value
      IF ( ANY(Frequency < MIN_FREQUENCY) ) THEN
        WRITE( msg,'("Input frequencies < ",f6.1," GHz found.")' ) MIN_FREQUENCY
        CALL Display_Message( ROUTINE_NAME, msg, WARNING )
      END IF
      IF ( ANY(Frequency > MAX_FREQUENCY) ) THEN
        WRITE( msg,'("Input frequencies > ",f6.1," GHz found.")' ) MAX_FREQUENCY
        CALL Display_Message( ROUTINE_NAME, msg, WARNING )
      END IF

      ! Temperatures < min value or > max_value
      IF ( Temperature < MIN_TEMPERATURE ) THEN
        WRITE( msg,'("Input temperature < ",f6.1," K.")' ) MIN_TEMPERATURE
        CALL Display_Message( ROUTINE_NAME, msg, WARNING )
      END IF
      IF ( Temperature > MAX_TEMPERATURE ) THEN
        WRITE( msg,'("Input temperature > ",f6.1," K.")' ) MAX_TEMPERATURE
        CALL Display_Message( ROUTINE_NAME, msg, WARNING )
      END IF

    END IF

    ! Calculate the relative reciprocal temperature, Theta,
    ! and common line strength factor, Common_s1
    Theta     = Compute_Theta( Temperature )
    Common_s1 = Compute_Common_s1( Theta )

    ! Calculate the wet line attentuation
    WetLine_Compute: IF ( Compute_WetLine ) THEN
      ! ...Calculate the H2O line strength factor
      H2O_s1 = Compute_H2O_s1( H2O_Pressure, Theta )
      ! ...Loop over the required frequencies
      WetLine_Frequency_Loop: DO l = 1, n_Frequencies
        ! ...Compute the water vapour line absoprtion
        H2O_Line_Absorption = Compute_H2O_Line_Absorption( &
                                Frequency(l)    , &
                                Dry_Air_Pressure, &
                                H2O_Pressure    , &
                                Theta           , &
                                Common_s1       , &
                                H2O_s1            )
        IF ( H2O_Line_Absorption < ZERO ) THEN
          IF ( noisy ) THEN
            WRITE( msg,'("Negative sum for H2O line absorption, ",es13.6,1x,&
                        &"at Frequency ",f8.3,1x,&
                        &"GHz. Setting to 0.0.")' ) H2O_Line_Absorption, Frequency(l)
            CALL Display_Message( ROUTINE_NAME, msg, WARNING )
          END IF
          H2O_Line_Absorption = ZERO
        END IF
        ! ...Compute the attenuation in dB/km
        WetLine_Attenuation(l) = ATT_K * Frequency(l) * H2O_Line_Absorption
      END DO WetLine_Frequency_Loop
    END IF WetLine_Compute

    ! Calculate the wet continuum attentuation for each frequency
    WetContinuum_Compute: IF ( Compute_WetContinuum ) THEN
      ! ...Loop over the required frequencies
      WetContinuum_Frequency_Loop: DO l = 1, n_Frequencies
!***TURNING OFF WET CONTINUUM***
!        ! ...Compute the water vapour continuum absoprtion
!        WetContinuum_Absorption = Compute_Wet_Continuum( &
!                                    Frequency(l)    , &
!                                    Dry_Air_Pressure, &
!                                    H2O_Pressure    , &
!                                    Theta             )
!        IF ( WetContinuum_Absorption < ZERO ) THEN
!          IF ( noisy ) THEN
!            WRITE( msg,'("Negative value for wet continuum absorption, ",es13.6,1x,&
!                        &"at Frequency ",f8.3,1x,&
!                        &"GHz. Setting to 0.0.")' ) WetContinuum_Absorption, Frequency(l)
!            CALL Display_Message( ROUTINE_NAME, msg, WARNING )
!          END IF
!          WetContinuum_Absorption = ZERO
!        END IF
WetContinuum_Absorption = ZERO
!***TURNING OFF WET CONTINUUM***
        ! ...Compute the attenuation in dB/km
        WetContinuum_Attenuation(l) = ATT_K * Frequency(l) * WetContinuum_Absorption
      END DO WetContinuum_Frequency_Loop
    END IF WetContinuum_Compute

    ! Calculate the dry line attentuation for each frequency
    DryLine_Compute: IF ( Compute_DryLine ) THEN
      ! ...Calculate the O2 line strength, width, and overlap factors
      O2_s1 = Compute_O2_s1( Dry_Air_Pressure, Theta )
      O2_g1 = Compute_O2_g1( H2O_Pressure, Theta )
      O2_d1 = Compute_O2_d1( Dry_Air_Pressure + H2O_Pressure, Theta )
      ! ...Loop over the required frequencies
      DryLine_Frequency_Loop: DO l = 1, n_Frequencies
        ! ...Compute the oxygen line absoprtion
        O2_Line_Absorption = Compute_O2_Line_Absorption( &
                               Frequency(l)    , &
                               Dry_Air_Pressure, &
                               Theta           , &
                               Common_s1       , &
                               O2_s1           , &
                               O2_g1           , &
                               O2_d1             )
        IF ( O2_Line_Absorption < ZERO ) THEN
          IF ( noisy ) THEN
            WRITE( msg,'("Negative sum for O2 line absorption, ",es13.6,1x,&
                        &" at Frequency ",f8.3,1x,&
                        &"GHz. Setting to 0.0.")' ) O2_Line_Absorption, Frequency(l)
            CALL Display_Message( ROUTINE_NAME, msg, WARNING )
          END IF
          O2_Line_Absorption = ZERO
        END IF
        ! ...Compute the attenuation in dB/km
        DryLine_Attenuation(l) = ATT_K * Frequency(l) * O2_Line_Absorption
      END DO DryLine_Frequency_Loop
    END IF DryLine_Compute

    ! Calculate the dry continuum attentuation for each frequency
    DryContinuum_Compute: IF ( Compute_DryContinuum ) THEN
      ! ...Loop over the required frequencies
      DryContinuum_Frequency_Loop: DO l = 1, n_Frequencies
        ! ...Compute the non-resonant dry air absorption
        DryContinuum_Absorption = Compute_Dry_Continuum( &
                                    Frequency(l)    , &
                                    Dry_Air_Pressure, &
                                    H2O_Pressure    , &
                                    Theta             )
        IF ( DryContinuum_Absorption < ZERO ) THEN
          IF ( noisy ) THEN
            WRITE( msg,'("Negative sum for dry continuum absorption, ",es13.6,1x,&
                        &" at Frequency ",f8.3,1x,&
                        &"GHz. Setting to 0.0.")' ) DryContinuum_Absorption, Frequency(l)
            CALL Display_Message( ROUTINE_NAME, msg, WARNING )
          END IF
          DryContinuum_Absorption = ZERO
        END IF
        ! ...Compute the attenuation in dB/km
        DryContinuum_Attenuation(l) = ATT_K * Frequency(l) * DryContinuum_Absorption
      END DO DryContinuum_Frequency_Loop
    END IF DryContinuum_Compute

  END FUNCTION Liebe93_By_Frequency

  ! --------------------------------------
  ! Single-Frequency, Multi-Layer function
  ! --------------------------------------
  FUNCTION Liebe93_By_Layer( &
    Frequency               , &  ! Input
    Dry_Air_Pressure        , &  ! Input
    H2O_Pressure            , &  ! Input
    Temperature             , &  ! Input
    Quiet,                    &  ! Optional input
    WetLine_Attenuation     , &  ! Optional output
    WetContinuum_Attenuation, &  ! Optional output
    DryLine_Attenuation     , &  ! Optional output
    DryContinuum_Attenuation) &  ! Optional output
  RESULT( err_stat )
    ! Arguments
    REAL(fp),           INTENT(IN)  :: Frequency
    REAL(fp),           INTENT(IN)  :: Dry_Air_Pressure(:)
    REAL(fp),           INTENT(IN)  :: H2O_Pressure(:)
    REAL(fp),           INTENT(IN)  :: Temperature(:)
    LOGICAL , OPTIONAL, INTENT(IN)  :: Quiet
    REAL(fp), OPTIONAL, INTENT(OUT) :: WetLine_Attenuation(:)
    REAL(fp), OPTIONAL, INTENT(OUT) :: WetContinuum_Attenuation(:)
    REAL(fp), OPTIONAL, INTENT(OUT) :: DryLine_Attenuation(:)
    REAL(fp), OPTIONAL, INTENT(OUT) :: DryContinuum_Attenuation(:)
    ! Function result
    INTEGER :: err_stat
    ! Local parameters
    CHARACTER(*),  PARAMETER :: ROUTINE_NAME = 'Liebe93::By_Layer'
    ! Local variables
    CHARACTER(ML) :: msg
    LOGICAL :: noisy
    LOGICAL :: Compute_WetLine
    LOGICAL :: Compute_WetContinuum
    LOGICAL :: Compute_DryLine
    LOGICAL :: Compute_DryContinuum
    INTEGER :: n_Layers, k
    REAL(fp) :: Theta
    REAL(fp) :: O2_s1, Common_s1, O2_g1, O2_d1, H2O_s1
    REAL(fp) :: H2O_Line_Absorption
    REAL(fp) :: WetContinuum_Absorption
    REAL(fp) :: O2_Line_Absorption
    REAL(fp) :: DryContinuum_Absorption

    ! Setup
    err_stat = SUCCESS
    ! ...Check for consistent input array sizes
    n_Layers = SIZE(Dry_Air_Pressure)
    IF ( SIZE(H2O_Pressure) /= n_Layers .OR. &
         SIZE(Temperature ) /= n_Layers      ) THEN
      err_stat = FAILURE
      msg = 'Input argument arrays have inconsistent sizes.'
      CALL Display_Message( ROUTINE_NAME, msg, err_stat ); RETURN
    END IF
    ! ...Process keywords
    noisy = .TRUE.
    IF ( PRESENT(Quiet) ) noisy = .NOT.Quiet

    ! Check output argument validity
    ! ...The water vapour line attenuation
    Compute_WetLine = .FALSE.
    IF ( PRESENT(WetLine_Attenuation) ) THEN
      Compute_WetLine = .TRUE.
      IF ( SIZE(WetLine_Attenuation) /= n_layers ) THEN
        err_stat = FAILURE
        msg = 'WetLine_Attenuation argument array has inconsistent size.'
        CALL Display_Message( ROUTINE_NAME, msg, err_stat ); RETURN
      END IF
    END IF
    ! ...The water vapour continuum attenuation
    Compute_WetContinuum = .FALSE.
    IF ( PRESENT(WetContinuum_Attenuation) ) THEN
      Compute_WetContinuum = .TRUE.
      IF ( SIZE(WetContinuum_Attenuation) /= n_layers ) THEN
        err_stat = FAILURE
        msg = 'WetContinuum_Attenuation argument array has inconsistent size.'
        CALL Display_Message( ROUTINE_NAME, msg, err_stat ); RETURN
      END IF
    END IF
    ! ...The water vapour line attenuation
    Compute_DryLine = .FALSE.
    IF ( PRESENT(DryLine_Attenuation) ) THEN
      Compute_DryLine = .TRUE.
      IF ( SIZE(DryLine_Attenuation) /= n_layers ) THEN
        err_stat = FAILURE
        msg = 'DryLine_Attenuation argument array has inconsistent size.'
        CALL Display_Message( ROUTINE_NAME, msg, err_stat ); RETURN
      END IF
    END IF
    ! ...The water vapour continuum attenuation
    Compute_DryContinuum = .FALSE.
    IF ( PRESENT( DryContinuum_Attenuation ) ) THEN
      Compute_DryContinuum = .TRUE.
      IF ( SIZE( DryContinuum_Attenuation ) /= n_layers ) THEN
        err_stat = FAILURE
        msg = 'DryContinuum_Attenuation argument array has inconsistent size.'
        CALL Display_Message( ROUTINE_NAME, msg, err_stat ); RETURN
      END IF
    END IF
    ! ...If no output arguments at all, return
    IF ( ( .NOT. Compute_WetLine      ) .AND. &
         ( .NOT. Compute_WetContinuum ) .AND. &
         ( .NOT. Compute_DryLine      ) .AND. &
         ( .NOT. Compute_DryContinuum )       ) RETURN
    ! ...Check for invalid frequency
    IF ( Frequency < ZERO ) THEN
      err_stat = FAILURE
        msg = 'Input frequency must be > 0.0.'
        CALL Display_Message( ROUTINE_NAME, msg, err_stat ); RETURN
    END IF
    ! ...Check for invalid profile data
    IF ( ANY(Dry_Air_Pressure < TOLERANCE) .OR. &
         ANY(H2O_Pressure     < TOLERANCE) .OR. &
         ANY(Temperature      < TOLERANCE)      ) THEN
      err_stat = FAILURE
      msg = 'Input tressure/temperature arguments must be > 0.0.'
      CALL Display_Message( ROUTINE_NAME, msg, err_stat ); RETURN
    END IF

    ! Issue warnings for data outside range
    IF ( noisy )  THEN

      ! Frequency < min value or > max_value
      IF ( Frequency < MIN_FREQUENCY ) THEN
        WRITE( msg,'("Input frequency < ",f6.1," GHz.")' ) MIN_FREQUENCY
        CALL Display_Message( ROUTINE_NAME, msg, WARNING )
      END IF
      IF ( Frequency > MAX_FREQUENCY ) THEN
        WRITE( msg,'("Input frequency > ",f6.1," GHz.")' ) MAX_FREQUENCY
        CALL Display_Message( ROUTINE_NAME, msg, WARNING )
      END IF

      ! Temperatures < min value or > max_value
      IF ( ANY(Temperature < MIN_TEMPERATURE) ) THEN
        WRITE( msg,'("Input temperatures < ",f6.1," K found.")' ) MIN_TEMPERATURE
        CALL Display_Message( ROUTINE_NAME, msg, WARNING )
      END IF
      IF ( ANY(Temperature > MAX_TEMPERATURE) ) THEN
        WRITE( msg,'("Input temperatures > ",f6.1," K found.")' ) MAX_TEMPERATURE
        CALL Display_Message( ROUTINE_NAME, msg, WARNING )
      END IF

    END IF

    ! Begin the loop over atmospheric layers
    Layer_Loop: DO k = 1, n_Layers

      ! Calculate the reciprocal temperature, theta,
      ! and common line strength factor, Common_s1
      Theta     = Compute_Theta( Temperature(k) )
      Common_s1 = Compute_Common_s1( Theta )

      ! Calculate the wet line attentuation
      WetLine_Compute: IF ( Compute_WetLine ) THEN
        ! ...Calculate the H2O absorption equation coefficients
        H2O_s1 = Compute_H2O_s1( H2O_Pressure(k), Theta )
        ! ...Water vapour line absorption
        H2O_Line_Absorption = Compute_H2O_line_Absorption( &
                                Frequency          , &
                                Dry_Air_Pressure(k), &
                                H2O_Pressure(k)    , &
                                Theta              , &
                                Common_s1          , &
                                H2O_s1               )
        IF ( H2O_Line_Absorption < ZERO ) THEN
          IF ( noisy ) THEN
            WRITE( msg,'("Negative sum for H2O line absorption, ",es13.6,1x,&
                        &"in layer ",i0," for frequency ",f8.3,1x, &
                        &"GHz. Setting to 0.0.")' ) H2O_Line_Absorption, k, Frequency
            CALL Display_Message( ROUTINE_NAME, msg, WARNING )
          END IF
          H2O_Line_Absorption = ZERO
        END IF
        ! ...Calculate the attenuation in dB/km
        WetLine_Attenuation(k) = ATT_K * Frequency * H2O_Line_Absorption
      END IF WetLine_Compute

      ! Calculate the wet continuum attentuation
      WetContinuum_Compute: IF ( Compute_WetContinuum ) THEN
!***TURNING OFF WET CONTINUUM***
!        ! ...Compute the water vapour continuum
!        WetContinuum_Absorption = Compute_Wet_Continuum( &
!                                    Frequency          , &
!                                    Dry_Air_Pressure(k), &
!                                    H2O_Pressure(k)    , &
!                                    Theta                )
!        IF ( WetContinuum_Absorption < ZERO ) THEN
!          IF ( noisy ) THEN
!            WRITE( msg,'("Negative value for wet continuum absorption, ",es13.6,1x,&
!                        &"in layer ",i0," for frequency ",f8.3,1x, &
!                        &"GHz. Setting to 0.0.")' ) WetContinuum_Absorption, k, Frequency
!            CALL Display_Message( ROUTINE_NAME, msg, WARNING )
!          END IF
!          WetContinuum_Absorption = ZERO
!        END IF
WetContinuum_Absorption = ZERO
!***TURNING OFF WET CONTINUUM***
        ! ...Compute the attenuation in dB/km
        WetContinuum_Attenuation(k) = ATT_K * Frequency * WetContinuum_Absorption
      END IF WetContinuum_Compute

      ! Calculate the dry line attentuation
      DryLine_Compute: IF ( Compute_DryLine ) THEN
        ! ...Calculate the O2 line strength, width, and overlap factors
        O2_s1 = Compute_O2_s1( Dry_Air_Pressure(k), Theta )
        O2_g1 = Compute_O2_g1( H2O_Pressure(k), Theta )
        O2_d1 = Compute_O2_d1( Dry_Air_Pressure(k) + H2O_Pressure(k), Theta )
        ! ...Summation of the absorption lines of oxygen
        O2_Line_Absorption = Compute_O2_Line_Absorption( &
                               Frequency          , &
                               Dry_Air_Pressure(k), &
                               Theta              , &
                               Common_s1          , &
                               O2_s1              , &
                               O2_g1              , &
                               O2_d1                )
        IF ( O2_Line_Absorption < ZERO ) THEN
          IF ( noisy ) THEN
            WRITE( msg,'("Negative value for O2 line absorption, ",es13.6,1x,&
                        &"in layer ",i0," for frequency ",f8.3,1x, &
                        &"GHz. Setting to 0.0.")' ) O2_Line_Absorption, k, Frequency
            CALL Display_Message( ROUTINE_NAME, msg, WARNING )
          END IF
          O2_Line_Absorption = ZERO
        END IF
        ! ...Compute the attenuation in dB/km
        DryLine_Attenuation(k) = ATT_K * Frequency * O2_Line_Absorption
      END IF DryLine_Compute

      ! Calculate the dry continuum attentuation
      DryContinuum_Compute: IF ( Compute_DryContinuum ) THEN
        ! ...Compute the non-resonant dry air absorption
        !    Formula from Liebe 89 paper + correction of Liebe himself
        DryContinuum_Absorption = Compute_Dry_Continuum( Frequency, &
                                                         Dry_Air_Pressure(k), &
                                                         H2O_Pressure(k), &
                                                         Theta )
        IF ( DryContinuum_Absorption < ZERO ) THEN
          IF ( noisy ) THEN
            WRITE( msg,'("Negative value for dry continuum absorption, ",es13.6,1x,&
                        &"in layer ",i0," for frequency ",f8.3,1x, &
                        &"GHz. Setting to 0.0.")' ) DryContinuum_Absorption, k, Frequency
            CALL Display_Message( ROUTINE_NAME, msg, WARNING )
          END IF
          DryContinuum_Absorption = ZERO
        END IF
        ! ...Compute the attenuation in dB/km
        DryContinuum_Attenuation(k) = ATT_K * Frequency * DryContinuum_Absorption
      END IF DryContinuum_Compute

    END DO Layer_Loop

  END FUNCTION Liebe93_By_Layer

!##################################################################################
!##################################################################################
!##                                                                              ##
!##                          ## PRIVATE MODULE ROUTINES ##                       ##
!##                                                                              ##
!##################################################################################
!##################################################################################

  ! --------------------------------
  ! Reference reciprocal temperature
  ! --------------------------------
  FUNCTION Compute_Theta( Temperature ) RESULT( Theta )
    REAL(fp), INTENT(IN) :: Temperature
    REAL(fp)             :: Theta
    REAL(fp), PARAMETER :: REFERENCE_TEMPERATURE = 300.0_fp
    Theta = REFERENCE_TEMPERATURE / Temperature
  END FUNCTION Compute_Theta

  ! --------------------------------------
  ! Common O2 and H2O line strength factor
  ! --------------------------------------
  FUNCTION Compute_Common_s1( Theta ) RESULT( Common_s1 )
    REAL(fp), INTENT(IN) :: Theta
    REAL(fp)             :: Common_s1
    Common_s1 = ONE - Theta
  END FUNCTION Compute_Common_s1

  ! -----------------------
  ! O2 line strength factor
  ! -----------------------
  FUNCTION Compute_O2_s1( Dry_Air_P, Theta ) RESULT( O2_s1 )
    REAL(fp), INTENT(IN) :: Dry_Air_P
    REAL(fp), INTENT(IN) :: Theta
    REAL(fp)             :: O2_s1
    O2_s1 = Dry_Air_P * ( Theta**THREE )
  END FUNCTION Compute_O2_s1

  ! --------------------
  ! O2 line width factor
  ! --------------------
  FUNCTION Compute_O2_g1( H2O_P, Theta ) RESULT( O2_g1 )
    REAL(fp), INTENT(IN) :: H2O_P
    REAL(fp), INTENT(IN) :: Theta
    REAL(fp)             :: O2_g1
    O2_g1 = ONEpointONE * H2O_P * Theta
  END FUNCTION Compute_O2_g1

  ! ----------------------
  ! O2 line overlap factor
  ! ----------------------
  FUNCTION Compute_O2_d1( Total_P, Theta ) RESULT( O2_d1 )
    REAL(fp), INTENT(IN) :: Total_P
    REAL(fp), INTENT(IN) :: Theta
    REAL(fp)             :: O2_d1
    O2_d1 = Total_P * ( Theta**ZEROpointEIGHT )
  END FUNCTION Compute_O2_d1

  ! ------------------------
  ! H2O line strength factor
  ! ------------------------
  FUNCTION Compute_H2O_s1( H2O_P, Theta ) RESULT( H2O_s1 )
    REAL(fp), INTENT(IN) :: H2O_P
    REAL(fp), INTENT(IN) :: Theta
    REAL(fp)             :: H2O_s1
    H2O_s1 = H2O_P * ( Theta**THREEpointFIVE )
  END FUNCTION Compute_H2O_s1

!------------------------------------------------------------------------------
!
! NAME:
!       Compute_O2_Line_Absorption
!
! PURPOSE:
!       Function to calculate the oxygen line terms of the imaginary part of
!       the atmospheric refractivity.
!
! CALLING SEQUENCE:
!       O2_Line_Absorption = Compute_O2_Line_Absorption( f,         &  ! Input
!                                                        pd,        &  ! Input
!                                                        Theta,     &  ! Input
!                                                        Common_s1, &  ! Input
!                                                        O2_s1,     &  ! Input
!                                                        O2_g1,     &  ! Input
!                                                        O2_d1      )  ! Input
!
! INPUT ARGUMENTS:
!       f:                   Frequency
!                            UNITS:      GHz
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Scalar
!                            ATTRIBUTES: INTENT(IN)
!
!       pd:                  Dry air partial pressure.
!                            UNITS:      hectoPascals, hPa
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Scalar
!                            ATTRIBUTES: INTENT(IN)
!
!       Theta:               Reciprocal temperature ratio, 300/T
!                            UNITS:      None.
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Scalar
!                            ATTRIBUTES: INTENT(IN)
!
!       Common_s1:           Line strength factor common to both O2 and H2O
!                            forumlations, (1 - Theta)
!                            UNITS:      None.
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Scalar
!                            ATTRIBUTES: INTENT(IN)
!
!       O2_s1:               O2 line strength factor, pd * Theta^3
!                            UNITS:      kPa
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Scalar
!                            ATTRIBUTES: INTENT(IN)
!
!       O2_g1:               O2 line width factor, 1.1 * pw * Theta
!                            UNITS:      kPa
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Scalar
!                            ATTRIBUTES: INTENT(IN)
!
!       O2_d1:               O2 line overlap factor, pd * Theta^0.8
!                            UNITS:      kPa
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Scalar
!                            ATTRIBUTES: INTENT(IN)
!
! FUNCTION RESULT:
!       O2_Line_Absorption:  Dry air refractivity due to O2 line
!                            absorption
!                            UNITS:      ppm
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Scalar
!
! PROCEDURE:
!       Liebe, H.J., G.A. Hufford, and M.G. Cotton, 1993: "Propagation
!         modelling of moist air and suspended water/ice particles at
!         frequencies below 1000 GHz". In AGARD Conference Proceedings
!         542, Atmospheric propagation effects through natural and man-made
!         obscurants for visible through MM-wave radiation, pp 3.1-3.10.
!         Presented at the Electromagnetic Wave Propagation Panel Symposium,
!         Palma de Mallorca, Spain, 17-20 May 1993.
!
!       Synopsis:
!
!       The imaginary part of the refractivity due to O2 line absorption in
!       moist air can be expressed by,
!               __
!              \
!         ND =  > S(f).F(f)    ppm  .....(2)
!              /__
!                i
!
!       with F being the line shape function derived from the complex function
!       shown in eqn.(3) of the MPM93 paper,
!
!                f   [    g - d.(fo - f)         g - d.(fo + f)    ]
!        F(f) = ---- [ -------------------- + -------------------- ]
!                fo  [  ( fo - f )^2 + g^2     ( fo + f )^2 + g^2  ]
!
!       where i = O2 line index
!             f = frequency
!             S = O2 line strength
!             g = O2 line width (gamma)
!             d = O2 line overlap (delta)
!
!       The line strength, width, and overlap are modeled using,
!
!         S = a1 . pd . theta^3 . EXP[ a2 . ( 1 - theta ) ]
!
!         g = a3.10^-3 ( pd.theta^a4 + 1.1.pw.theta )
!
!         d = ( a5 + a6.theta ) . pd . theta^0.8
!
!       with a1,a2 = line strength coefficients
!            a3,a4 = line width coefficients
!            a5,a6 = line overlap coefficients
!            pd    = dry gas partial pressure
!            pw    = water vapor partial pressure
!            theta = reciprocal temperature ratio, 300/T
!
!       A rough estimate of line behaviour in the mesosphere is obtained by
!       g with gh where
!                 _______________
!                /
!         gh = \/ g^2 + 625.B^2
!
!       where B = the magnetic field strength (22-65uT) depending on the
!                 geographic location and altitude
!
!       The Zeeman parameter approximates this solution.
!
!       All equation numbers refer to the referenced Liebe93 paper.
!
!
! CREATION HISTORY:
!       F77 version written by : L. Phalippou 13-Dec-1993
!
!       F90 version translation: Paul van Delst, CIMSS/SSEC 21-Apr-2002
!                                paul.vandelst@ssec.wisc.edu
!
!------------------------------------------------------------------------------

  FUNCTION Compute_O2_Line_Absorption( f,         &  ! Input
                                       pd,        &  ! Input
                                       Theta,     &  ! Input
                                       Common_s1, &  ! Input
                                       O2_s1,     &  ! Input
                                       O2_g1,     &  ! Input
                                       O2_d1 )    &  ! Input
                                     RESULT( O2_Line_Absorption )
    ! Arguments
    REAL(fp), INTENT(IN) :: f
    REAL(fp), INTENT(IN) :: pd
    REAL(fp), INTENT(IN) :: Theta
    REAL(fp), INTENT(IN) :: Common_s1
    REAL(fp), INTENT(IN) :: O2_s1
    REAL(fp), INTENT(IN) :: O2_g1
    REAL(fp), INTENT(IN) :: O2_d1
    ! Function result
    REAL(fp) :: O2_Line_Absorption
    ! Local parameters
    REAL(fp), PARAMETER :: ZEEMAN_K = ( 25.0_fp * 0.6e-04_fp )**2
    ! Local variables
    INTEGER :: i
    REAL(fp) :: Si
    REAL(fp) :: Gamma, Gamma2
    REAL(fp) :: Delta
    REAL(fp) :: f_diff, f_sum
    REAL(fp) :: Fi_n1, Fi_n2
    REAL(fp) :: Fi_d1, Fi_d2
    REAL(fp) :: Fi

    ! Initialise the return value
    ! ---------------------------
    O2_Line_Absorption = ZERO

    ! Loop over the oxygen lines
    ! --------------------------
    DO i = 1, N_O2_LINES

      ! Calculate the line strength
      ! ---------------------------
      Si = O2_A1(i) * O2_s1 * EXP( O2_A2(i) * Common_s1 )

      ! Calculate the line width and overlap parameters
      ! -----------------------------------------------
      ! Calculate the pressure-broadened line width
      Gamma  = O2_A3(i) * ( ( pd * ( Theta**(ZEROpointEIGHT - O2_A4(i)) ) ) + O2_g1 )
      ! Zeeman effect in the mesosphere due to magnetic field
!!      Gamma  = SQRT( Gamma**2 + ZEEMAN_K )
      ! Compute the square of the line width
      Gamma2 = Gamma * Gamma
      ! Calculate the line overlap (interference)
      Delta = ( O2_A5(i) + ( O2_A6(i) * Theta ) ) * O2_d1

      ! Calculate the line shape, F
      ! ---------------------------
      ! Frequency difference and sum from the line frequency
      f_diff = O2_LINE_FREQUENCY(i) - f
      f_sum  = O2_LINE_FREQUENCY(i) + f
      ! The numerator of the line shape function components
      Fi_n1 = Gamma - ( Delta * f_diff )
      Fi_n2 = Gamma - ( Delta * f_sum  )
      ! The denominator of the line shape components
      Fi_d1 = ( f_diff * f_diff ) + Gamma2
      Fi_d2 = ( f_sum  * f_sum  ) + Gamma2
      ! Calculate the absorption line shape
      Fi = ( f / O2_LINE_FREQUENCY(i) ) * ( ( Fi_n1/Fi_d1 ) + ( Fi_n2/Fi_d2 ) )

      ! Calculate the imaginary component of the
      ! dry air refractivity - the absorption
      ! ----------------------------------------
      O2_Line_Absorption = O2_Line_Absorption + ( Si * Fi )

    END DO

  END FUNCTION Compute_O2_Line_Absorption

!------------------------------------------------------------------------------
!
! NAME:
!       Compute_Dry_Continuum
!
! PURPOSE:
!       Function to calculate the nonresonant imaginary term of the dry air
!       refractivity.
!
! CALLING SEQUENCE:
!       Dry_Continuum = Compute_Dry_Continuum( f,         &  ! Input
!                                              pd,        &  ! Input
!                                              pw,        &  ! Input
!                                              Theta )    &  ! Input
!
! INPUT ARGUMENTS:
!       f:                 Frequency
!                          UNITS:      GHz
!                          TYPE:       REAL(fp)
!                          DIMENSION:  Scalar
!                          ATTRIBUTES: INTENT(IN)
!
!       pd:                Dry air partial pressure.
!                          UNITS:      hectoPascals, hPa
!                          TYPE:       REAL(fp)
!                          DIMENSION:  Scalar
!                          ATTRIBUTES: INTENT(IN)
!
!       pw:                Water vapor partial pressure.
!                          UNITS:      hectoPascals, hPa
!                          TYPE:       REAL(fp)
!                          DIMENSION:  Scalar
!                          ATTRIBUTES: INTENT(IN)
!
!       Theta:             Reciprocal temperature ratio, 300/T
!                          UNITS:      None.
!                          TYPE:       REAL(fp)
!                          DIMENSION:  Scalar
!                          ATTRIBUTES: INTENT(IN)
!
! FUNCTION RESULT:
!       Dry_Continuum:     Nonresonant refractivity term for dry air
!                          UNITS:      ppm
!                          TYPE:       REAL(fp)
!                          DIMENSION:  Scalar
!
! PROCEDURE:
!       Liebe, H.J., G.A. Hufford, and M.G. Cotton, 1993: "Propagation
!         modelling of moist air and suspended water/ice particles at
!         frequencies below 1000 GHz". In AGARD Conference Proceedings
!         542, Atmospheric propagation effects through natural and man-made
!         obscurants for visible through MM-wave radiation, pp 3.1-3.10.
!         Presented at the Electromagnetic Wave Propagation Panel Symposium,
!         Palma de Mallorca, Spain, 17-20 May 1993.
!
!       Synopsis:
!
!       Nonresonant refractivity of dry air makes a small contribution at
!       surface pressures due to the Debye spectrum of oxygen below 10GHz and
!       the pressure-induced nitrogen absorption that becomes effective above
!       100GHz. The functional form is,
!
!         Nn = So.Fo(f) + i.Sn.Fn"(f)     ppm
!
!       where f     = frequency,
!             So,Fo = nonresonant O2 spectrum contribution,
!             Sn,Fn = pressure-induced N2 absorption contribution
!
!       and the imaginary part of Nn provides the attenuated power due to the
!       dry continuum term.
!
!
!       The non-resonant O2 strength and lineshape are given by,
!
!         So = 6.14x10^-5 . pd . theta^2
!
!       and
!
!                  -f
!         Fo = ----------
!               f + i.go
!
!       where go = relaxation frequency
!                = 0.56x10^-3.p.theta^0.8
!             p  = Total pressure
!
!       The imaginary part of the non-resonant lineshape is therefore,
!
!                  f.go
!         Fo = ------------
!               f^2 + go^2
!
!
!       The strength and line shape of the N2 absorption term is given by,
!
!         Sn = 1.40x10^-12 . pd^2 . theta^3.5
!
!       and
!
!                       f
!         Fn ~ ---------------------------
!               1 + ( 1.9x10^-5 . f^1.5 )
!
!
! CREATION HISTORY:
!       F77 version written by : L. Phalippou 13-Dec-1993
!
!       F90 version translation: Paul van Delst, CIMSS/SSEC 21-Apr-2002
!                                paul.vandelst@ssec.wisc.edu
!
!
!------------------------------------------------------------------------------

  FUNCTION Compute_Dry_Continuum( f,       &  ! Input
                                  pd,      &  ! Input
                                  pw,      &  ! Input
                                  Theta )  &  ! Input
                                RESULT( Dry_Continuum )
    ! Arguments
    REAL(fp), INTENT(IN) :: f
    REAL(fp), INTENT(IN) :: pd
    REAL(fp), INTENT(IN) :: pw
    REAL(fp), INTENT(IN) :: Theta
    ! Function result
    REAL(fp) :: Dry_Continuum
    ! Parameters
    REAL(fp), PARAMETER :: So_K    = 6.14e-05_fp
    REAL(fp), PARAMETER :: GAMMA_K = 5.6e-04_fp
    REAL(fp), PARAMETER :: Sn_K    = 1.4e-12_fp
    REAL(fp), PARAMETER :: Fn_K    = 1.93e-05_fp
    ! Local variables
    REAL(fp) :: Gamma
    REAL(fp) :: So, Fo
    REAL(fp) :: Sn, Fn

    ! Calculate the effect of the non-resonant oxygen spectrum
    ! ...Calculate the strength
    So = So_K * pd * ( Theta**2 )
    ! ...Calculate the width
    Gamma = GAMMA_K * ( pd + pw ) * ( Theta**ZEROpointEIGHT )
    ! ...Calculate the lineshape
    Fo = ( f * Gamma ) / ( f**2 + Gamma**2 )

    ! Calculate the effect of pressure induced N2 absorption
    ! ...Calculate the N2 absorption strength
    Sn = Sn_K * ( pd**2 ) * ( theta**THREEpointFIVE )
    ! ...Calculate the shape
    Fn = f / ( ONE + ( Fn_K * ( f**ONEpointFIVE ) ) )

    ! Sum the components
    Dry_Continuum = ( So * Fo ) + ( Sn * Fn )

  END FUNCTION Compute_Dry_Continuum

!------------------------------------------------------------------------------
!
! NAME:
!       Compute_H2O_Line_Absorption
!
! PURPOSE:
!       Function to calculate the water vapor line terms of the imaginary part
!       of the atmospheric refractivity.
!
!       A continuum term is included by considering a pseudo-line centred out
!       of band at 1780GHz. The parameters for this pseudo-line are in the
!       coefficient data modules.
!
! CALLING SEQUENCE:
!       H2O_Line_Absorption = Compute_H2O_Line_Absorption( f,         &  ! Input
!                                                          pd,        &  ! Input
!                                                          pw,        &  ! Input
!                                                          Theta,     &  ! Input
!                                                          Common_s1, &  ! Input
!                                                          H2O_s1 )   &  ! Input
!
! INPUT ARGUMENTS:
!       f:                 Frequency
!                          UNITS:      GHz
!                          TYPE:       REAL(fp)
!                          DIMENSION:  Scalar
!                          ATTRIBUTES: INTENT(IN)
!
!       pd:                Dry air partial pressure.
!                          UNITS:      kiloPascals, kPa
!                          TYPE:       REAL(fp)
!                          DIMENSION:  Scalar
!                          ATTRIBUTES: INTENT(IN)
!
!       pw:                Water vapor partial pressure.
!                          UNITS:      kiloPascals, kPa
!                          TYPE:       REAL(fp)
!                          DIMENSION:  Scalar
!                          ATTRIBUTES: INTENT(IN)
!
!       Theta:             Reciprocal temperature ratio, 300/T
!                          UNITS:      None.
!                          TYPE:       REAL(fp)
!                          DIMENSION:  Scalar
!                          ATTRIBUTES: INTENT(IN)
!
!       Common_s1:         Line strength factor common to both H2O and H2O
!                          forumlations, (1 - Theta)
!                          UNITS:      None.
!                          TYPE:       REAL(fp)
!                          DIMENSION:  Scalar
!                          ATTRIBUTES: INTENT(IN)
!
!       H2O_s1:            H2O line strength factor, pw * Theta^3.5
!                          UNITS:      kPa
!                          TYPE:       REAL(fp)
!                          DIMENSION:  Scalar
!                          ATTRIBUTES: INTENT(IN)
!
! FUNCTION RESULT:
!       H2O_Line_Absorption: Atmospheric refractivity due to H2O line
!                            absorption
!                            UNITS:      ppm
!                            TYPE:       REAL(fp)
!                            DIMENSION:  Scalar
!
! PROCEDURE:
!       Liebe, H.J., G.A. Hufford, and M.G. Cotton, 1993: "Propagation
!         modelling of moist air and suspended water/ice particles at
!         frequencies below 1000 GHz". In AGARD Conference Proceedings
!         542, Atmospheric propagation effects through natural and man-made
!         obscurants for visible through MM-wave radiation, pp 3.1-3.10.
!         Presented at the Electromagnetic Wave Propagation Panel Symposium,
!         Palma de Mallorca, Spain, 17-20 May 1993.
!
!       Synopsis:
!
!       The imaginary part of the refractivity due to H2O line absorption in
!       moist air can be expressed by,
!               __
!              \
!         NW =  > S(f).F(f)    ppm  .....(5)
!              /__
!                i
!
!       with F being the line shape function derived from the complex function
!       shown in eqn.(3) of the MPM93 paper,
!
!                f   [    g - d.(fo - f)         g - d.(fo + f)    ]
!        F(f) = ---- [ -------------------- + -------------------- ]
!                fo  [  ( fo - f )^2 + g^2     ( fo + f )^2 + g^2  ]
!
!       where i = H2O line index
!             f = frequency
!             S = H2O line strength
!             g = H2O line width (gamma)
!             d = H2O line overlap (delta)
!
!       The line strength, width, and overlap are modeled using,
!
!         S = a1 . pw . theta^3.5 . EXP[ a2 . ( 1 - theta ) ]
!
!         g = a3.10^-3 ( a4.pw.theta^a6 + pd.theta^a5 )
!
!         d = 0
!
!       with a1,a2       = line strength coefficients
!            a3,a4,a5,a6 = line width coefficients
!            pd          = dry gas partial pressure
!            pw          = water vapor partial pressure
!            theta       = reciprocal temperature ratio, 300/T
!
!       Line overlap (delta) is neglected (d(i) = 0) and Doppler-broadening
!       is approximated for pressures below 0.7mb by
!                                  _______________________
!                                 /
!         g*(i) = 0.535.g(i) +   /  0.217.g(i)^2 + g(D)^2
!                              \/
!
!       where g(D) = Doppler width
!
!                     1.46x10^-6 . fo(i)
!                  = --------------------
!                          ______
!                         /
!                       \/ theta
!
!       Note that the water vapor continuum is modeled by considering a pseudo-
!       line centred out of band above 1 THz. Thus there is an extra line in the
!       coefficient data for a line centred at 1780GHz.
!
!
! CREATION HISTORY:
!       F77 version written by : L. Phalippou 13-Dec-1993
!
!       F90 version translation: Paul van Delst, CIMSS/SSEC 21-Apr-2002
!                                paul.vandelst@ssec.wisc.edu
!
!
!------------------------------------------------------------------------------

  FUNCTION Compute_H2O_Line_Absorption( f,         &  ! Input
                                        pd,        &  ! Input
                                        pw,        &  ! Input
                                        Theta,     &  ! Input
                                        Common_s1, &  ! Input
                                        H2O_s1 )   &  ! Input
                                      RESULT( H2O_Line_Absorption )
    ! Arguments
    REAL(fp), INTENT(IN) :: f
    REAL(fp), INTENT(IN) :: pd
    REAL(fp), INTENT(IN) :: pw
    REAL(fp), INTENT(IN) :: Theta
    REAL(fp), INTENT(IN) :: Common_s1
    REAL(fp), INTENT(IN) :: H2O_s1
    ! Function result
    REAL(fp) :: H2O_Line_Absorption
    ! Local parameters
    REAL(fp), PARAMETER :: GAMMA_D_K1 = 1.46e-06_fp
    REAL(fp), PARAMETER :: GAMMA_D_K2 = 0.535_fp
    REAL(fp), PARAMETER :: GAMMA_D_K3 = 0.217_fp
    ! Local variables
    INTEGER :: i
    REAL(fp) :: Si
    REAL(fp) :: Gamma, Gamma2, GammaD2
    REAL(fp) :: f_diff, f_sum
    REAL(fp) :: Fi_d1, Fi_d2
    REAL(fp) :: Fi

    ! Initialise the return value
    H2O_Line_Absorption = ZERO

    ! Loop over water vapour lines
    DO i = 1, N_H2O_LINES

      ! Calculate the line strength
      Si = H2O_A1(i) * H2O_s1 * EXP( H2O_A2(i) * Common_s1 )

      ! Calculate the line width
      ! ...Calculate the pressure-broadened width
      Gamma = H2O_A3(i) * ( (             pd * ( Theta**H2O_A5(i) ) ) + &
                            ( H2O_A4(i) * pw * ( Theta**H2O_A6(i) ) )   )
      ! ...Calculate the square of the Doppler width
      GammaD2 = ( ( GAMMA_D_K1 * H2O_LINE_FREQUENCY(i) )**2 ) / Theta
      ! ...Approximate the Doppler broadening
      Gamma   = ( GAMMA_D_K2 * Gamma ) + SQRT( ( GAMMA_D_K3 * Gamma**2 ) + GammaD2 )
      Gamma2  = Gamma * Gamma

      ! Calculate the line shape, F
      ! ...Frequency difference and sum from the line frequency
      f_diff = H2O_LINE_FREQUENCY(i) - f
      f_sum  = H2O_LINE_FREQUENCY(i) + f
      ! ...The denominator of the line shape components
      Fi_d1 = ( f_diff * f_diff ) + Gamma2
      Fi_d2 = ( f_sum  * f_sum  ) + Gamma2
      ! ...Calculate the absorption line shape
      Fi = ( f / H2O_LINE_FREQUENCY(i) ) * ( ( Gamma/Fi_d1 ) + ( Gamma/Fi_d2 )   )

      ! Calculate the imaginary component of the
      ! water vapor refractivity - the absorption
      H2O_Line_Absorption = H2O_Line_Absorption + ( Si * Fi )

    END DO

  END FUNCTION Compute_H2O_Line_Absorption

END MODULE MWLBL_Liebe93

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm