source: LMDZ6/trunk/libf/phylmd/cosp/gases.f90 @ 5450

! $Revision: 88 $, $Date: 2013-11-13 15:08:38 +0100 (mer. 13 nov. 2013) $
! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/quickbeam/gases.f90 $
  function gases(PRES_mb,T,RH,f)
  implicit none
  
! Purpose:
!   Compute 2-way gaseous attenuation through a volume in microwave
!
! Inputs:
!   [PRES_mb]   pressure (mb) (hPa)
!   [T]         temperature (K)
!   [RH]        relative humidity (%)
!   [f]         frequency (GHz), < 300 GHz
!
! Returns:
!   2-way gaseous attenuation (dB/km)
!
! Reference:
!   Uses method of Liebe (1985)
!
! Created:
!   12/09/05  John Haynes (haynes@atmos.colostate.edu)
! Modified:
!   01/31/06  Port from IDL to Fortran 90

  integer, parameter :: &
  nbands_o2 = 48 ,&
  nbands_h2o = 30
  real*8, intent(in) :: PRES_mb, T, RH, f
  real*8 :: gases, th, e, p, sumo, gm0, a0, ap, term1, term2, term3, &
            bf, be, term4, npp
  real*8, dimension(nbands_o2) :: v0, a1, a2, a3, a4, a5, a6
  real*8, dimension(nbands_h2o) :: v1, b1, b2, b3
  real*8 :: e_th,one_th,pth3,eth35,aux1,aux2,aux3,aux4
  real*8 :: gm,delt,x,y,gm2
  real*8 :: fpp_o2,fpp_h2o,s_o2,s_h2o
  integer :: i
  
! // table1 parameters  v0, a1, a2, a3, a4, a5, a6  
  data v0/49.4523790,49.9622570,50.4742380,50.9877480,51.5033500, &
  52.0214090,52.5423930,53.0669060,53.5957480,54.1299999,54.6711570, &
  55.2213650,55.7838000,56.2647770,56.3378700,56.9681000,57.6124810, &
  58.3238740,58.4465890,59.1642040,59.5909820,60.3060570,60.4347750, &
  61.1505580,61.8001520,62.4112120,62.4862530,62.9979740,63.5685150, &
  64.1277640,64.6789000,65.2240670,65.7647690,66.3020880,66.8368270, &
  67.3695950,67.9008620,68.4310010,68.9603060,69.4890210,70.0173420, &
  118.7503410,368.4983500,424.7631200,487.2493700,715.3931500, &
  773.8387300, 834.1453300/
  data a1/0.0000001,0.0000003,0.0000009,0.0000025,0.0000061,0.0000141, &
  0.0000310,0.0000641,0.0001247,0.0002280,0.0003918,0.0006316,0.0009535, &
  0.0005489,0.0013440,0.0017630,0.0000213,0.0000239,0.0000146,0.0000240, &
  0.0000211,0.0000212,0.0000246,0.0000250,0.0000230,0.0000193,0.0000152, &
  0.0000150,0.0000109,0.0007335,0.0004635,0.0002748,0.0001530,0.0000801, &
  0.0000395,0.0000183,0.0000080,0.0000033,0.0000013,0.0000005,0.0000002, &
  0.0000094,0.0000679,0.0006380,0.0002350,0.0000996,0.0006710,0.0001800/
  data a2/11.8300000,10.7200000,9.6900000,8.8900000,7.7400000,6.8400000, &
  6.0000000,5.2200000,4.4800000,3.8100000,3.1900000,2.6200000,2.1150000, &
  0.0100000,1.6550000,1.2550000,0.9100000,0.6210000,0.0790000,0.3860000, &
  0.2070000,0.2070000,0.3860000,0.6210000,0.9100000,1.2550000,0.0780000, &
  1.6600000,2.1100000,2.6200000,3.1900000,3.8100000,4.4800000,5.2200000, &
  6.0000000,6.8400000,7.7400000,8.6900000,9.6900000,10.7200000,11.8300000, &
  0.0000000,0.0200000,0.0110000,0.0110000,0.0890000,0.0790000,0.0790000/
  data a3/0.0083000,0.0085000,0.0086000,0.0087000,0.0089000,0.0092000, &
  0.0094000,0.0097000,0.0100000,0.0102000,0.0105000,0.0107900,0.0111000, &
  0.0164600,0.0114400,0.0118100,0.0122100,0.0126600,0.0144900,0.0131900, &
  0.0136000,0.0138200,0.0129700,0.0124800,0.0120700,0.0117100,0.0146800, &
  0.0113900,0.0110800,0.0107800,0.0105000,0.0102000,0.0100000,0.0097000, &
  0.0094000,0.0092000,0.0089000,0.0087000,0.0086000,0.0085000,0.0084000, &
  0.0159200,0.0192000,0.0191600,0.0192000,0.0181000,0.0181000,0.0181000/
  data a4/0.0000000,0.0000000,0.0000000,0.0000000,0.0000000,0.0000000, &
  0.0000000,0.0000000,0.0000000,0.0000000,0.0000000,0.0000000,0.0000000, &
  0.0000000,0.0000000,0.0000000,0.0000000,0.0000000,0.0000000,0.0000000, &
  0.0000000,0.0000000,0.0000000,0.0000000,0.0000000,0.0000000,0.0000000, &
  0.0000000,0.0000000,0.0000000,0.0000000,0.0000000,0.0000000,0.0000000, &
  0.0000000,0.0000000,0.0000000,0.0000000,0.0000000,0.0000000,0.0000000, &
  0.0000000,0.6000000,0.6000000,0.6000000,0.6000000,0.6000000,0.6000000/
  data a5/0.0056000,0.0056000,0.0056000,0.0055000,0.0056000,0.0055000, &
  0.0057000,0.0053000,0.0054000,0.0048000,0.0048000,0.0041700,0.0037500, &
  0.0077400,0.0029700,0.0021200,0.0009400,-0.0005500,0.0059700,-0.0024400, &
  0.0034400,-0.0041300,0.0013200,-0.0003600,-0.0015900,-0.0026600, &
  -0.0047700,-0.0033400,-0.0041700,-0.0044800,-0.0051000,-0.0051000, &
  -0.0057000,-0.0055000,-0.0059000,-0.0056000,-0.0058000,-0.0057000, &
  -0.0056000,-0.0056000,-0.0056000,-0.0004400,0.0000000,0.0000000, &
  0.0000000,0.0000000,0.0000000,0.0000000/
  data a6/1.7000000,1.7000000,1.7000000,1.7000000,1.8000000,1.8000000,&
  1.8000000,1.9000000,1.8000000,2.0000000,1.9000000,2.1000000,2.1000000, &
  0.9000000,2.3000000,2.5000000,3.7000000,-3.1000000,0.8000000,0.1000000, &
  0.5000000,0.7000000,-1.0000000,5.8000000,2.9000000,2.3000000,0.9000000, &
  2.2000000,2.0000000,2.0000000,1.8000000,1.9000000,1.8000000,1.8000000, &
  1.7000000,1.8000000,1.7000000,1.7000000,1.7000000,1.7000000,1.7000000, &
  0.9000000,1.0000000,1.0000000,1.0000000,1.0000000,1.0000000,1.0000000/

! // table2 parameters  v1, b1, b2, b3
  data v1/22.2350800,67.8139600,119.9959400,183.3101170,321.2256440, &
  325.1529190,336.1870000,380.1973720,390.1345080,437.3466670,439.1508120, &
  443.0182950,448.0010750,470.8889740,474.6891270,488.4911330,503.5685320, &
  504.4826920,556.9360020,620.7008070,658.0065000,752.0332270,841.0735950, &
  859.8650000,899.4070000,902.5550000,906.2055240,916.1715820,970.3150220, &
  987.9267640/
  data b1/0.1090000,0.0011000,0.0007000,2.3000000,0.0464000,1.5400000, &
  0.0010000,11.9000000,0.0044000,0.0637000,0.9210000,0.1940000,10.6000000, &
  0.3300000,1.2800000,0.2530000,0.0374000,0.0125000,510.0000000,5.0900000, &
  0.2740000,250.0000000,0.0130000,0.1330000,0.0550000,0.0380000,0.1830000, &
  8.5600000,9.1600000,138.0000000/
  data b2/2.1430000,8.7300000,8.3470000,0.6530000,6.1560000,1.5150000, &
  9.8020000,1.0180000,7.3180000,5.0150000,3.5610000,5.0150000,1.3700000, &
  3.5610000,2.3420000,2.8140000,6.6930000,6.6930000,0.1140000,2.1500000, &
  7.7670000,0.3360000,8.1130000,7.9890000,7.8450000,8.3600000,5.0390000, &
  1.3690000,1.8420000,0.1780000/
  data b3/0.0278400,0.0276000,0.0270000,0.0283500,0.0214000,0.0270000, &
  0.0265000,0.0276000,0.0190000,0.0137000,0.0164000,0.0144000,0.0238000, &
  0.0182000,0.0198000,0.0249000,0.0115000,0.0119000,0.0300000,0.0223000, &
  0.0300000,0.0286000,0.0141000,0.0286000,0.0286000,0.0264000,0.0234000, &
  0.0253000,0.0240000,0.0286000/
  
! // conversions
  th = 300./T       ! unitless
  e = (RH*th**5)/(41.45*10**(9.834*th-10))   ! kPa
  p = PRES_mb/10.-e ! kPa
  e_th = e*th
  one_th = 1 - th
  pth3 = p*th**(3)
  eth35 = e*th**(3.5)

! // term1
  sumo = 0.
  aux1 = 1.1*e_th
  do i=1,nbands_o2
    aux2 = f/v0(i)
    aux3 = v0(i)-f
    aux4 = v0(i)+f
    gm = a3(i)*(p*th**(0.8-a4(i))+aux1)
    gm2 = gm**2
    delt = a5(i)*p*th**a6(i)
    x = aux3**2+gm2
    y = aux4**2+gm2
    fpp_o2 = (((1./x)+(1./y))*(gm*aux2) - (delt*aux2)*((aux3/(x))-(aux4/(x))))
    s_o2 = a1(i)*pth3*exp(a2(i)*one_th)
    sumo = sumo + fpp_o2 * s_o2
  enddo
  term1 = sumo

! // term2
  gm0 = 5.6E-3*(p+1.1*e)*th**(0.8)
  a0 = 3.07E-4
  ap = 1.4*(1-1.2*f**(1.5)*1E-5)*1E-10
  term2 = (2*a0*(gm0*(1+(f/gm0)**2)*(1+(f/60.)**2))**(-1) + ap*p*th**(2.5)) &
          * f*p*th**2

! // term3
  sumo = 0.
  aux1 = 4.8*e_th
  do i=1,nbands_h2o
    aux2 = f/v1(i)
    aux3 = v1(i)-f
    aux4 = v1(i)+f
    gm = b3(i)*(p*th**(0.8)+aux1)
    gm2 = gm**2
    x = aux3**2+gm2
    y = aux4**2+gm2
    !delt = 0.
    fpp_h2o = ((1./x)+(1./y))*(gm*aux2) ! - (delt*aux2)*((aux3/(x))-(aux4/(x)))
    s_h2o = b1(i)*eth35*exp(b2(i)*one_th)
    sumo = sumo + fpp_h2o * s_h2o
  enddo
  term3 = sumo

! // term4
  bf = 1.4E-6
  be = 5.41E-5
  term4 = (bf*p+be*e*th**3)*f*e*th**(2.5)

! // summation and result
  npp = term1 + term2 + term3 + term4
  gases = 0.182*f*npp

  end function gases