!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE IRRAD                      ######
!######                                                      ######
!######                     Developed by                     ######
!######                                                      ######
!######    Goddard Cumulus Ensemble Modeling Group, NASA     ######
!######                                                      ######
!######     Center for Analysis and Prediction of Storms     ######
!######               University of Oklahoma                 ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE irrad(idim,jdim,m,n,np,                                      & 1,13
           taucl,ccld,pl,ta,wa,oa,co2,ts,                               &
           high,flx,flc,dfdts,st4,                                      &
           fclr,dbs,trant,th2o,tcon,tco2,                               &
           pa,dt,sh2o,swpre,swtem,sco3,scopre,scotem,                   &
           dh2o,dcont,dco2,do3,flxu,flxd,clr,blayer,                    &
           h2oexp,conexp,co2exp)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate IR fluxes due to water vapor, co2, and o3. Clouds in
!  different layers are assumed randomly overlapped.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: (a) Radiative Transfer Model: M.-D. Chou and M. Suarez
!          (b) Cloud Optics:Tao, Lang, Simpson, Sui, Ferrier and
!              Chou (1996)
!
!  MODIFICATION HISTORY:
!
!  03/15/1996 (Yuhe Liu)
!  Adopted the original code and formatted it in accordance with the
!  ARPS coding standard.
!
!-----------------------------------------------------------------------
!
!  ORIGINAL COMMENTS:
!
!******************** CLIRAD IR1  Date: Oct. 17, 1994 ****************
!*********************************************************************
!
! This routine computes ir fluxes due to water vapor, co2, and o3.
!   Clouds in different layers are assumed randomly overlapped.
!
! This is a vectorized code.  It computes fluxes simultaneously for
!   (m x n) soundings, which is a subset of (m x ndim) soundings.
!   In a global climate model, m and ndim correspond to the numbers of
!   grid boxes in the zonal and meridional directions, respectively.
!
! Detailed description of the radiation routine is given in
!   Chou and Suarez (1994).
!
! There are two options for computing cooling rate profiles.
!
!   if high = .true., transmission functions in the co2, o3, and the
!   three water vapor bands with strong absorption are computed using
!   table look-up.  cooling rates are computed accurately from the
!   surface up to 0.01 mb.
!   if high = .false., transmission functions are computed using the
!   k-distribution method with linear pressure scaling.  cooling rates
!   are not calculated accurately for pressures less than 20 mb.
!   the computation is faster with high=.false. than with high=.true.
!
! The IR spectrum is divided into eight bands:
!
!   bnad     wavenumber (/cm)   absorber         method
!
!    1           0 - 340           h2o            K/T
!    2         340 - 540           h2o            K/T
!    3         540 - 800       h2o,cont,co2       K,S,K/T
!    4         800 - 980       h2o,cont           K,S
!    5         980 - 1100      h2o,cont,o3        K,S,T
!    6        1100 - 1380      h2o,cont           K,S
!    7        1380 - 1900          h2o            K/T
!    8        1900 - 3000          h2o            K
!
! Note : "h2o" for h2o line absorption
!     "cont" for h2o continuum absorption
!     "K" for k-distribution method
!     "S" for one-parameter temperature scaling
!     "T" for table look-up
!
! The 15 micrometer region (540-800/cm) is further divided into
!   3 sub-bands :
!
!   subbnad   wavenumber (/cm)
!
!    1          540 - 620
!    2          620 - 720
!    3          720 - 800
!
!---- Input parameters                               units    size
!
!   number of soundings in zonal direction (m)        n/d      1
!   number of soundings in meridional direction (n)   n/d      1
!   maximum number of soundings in
!              meridional direction (ndim)         n/d      1
!   number of atmospheric layers (np)                 n/d      1
!   cloud optical thickness (taucl)                   n/d     m*ndim*np
!   cloud cover (ccld)                              fraction  m*ndim*np
!   level pressure (pl)                               mb      m*ndim*(np+1)
!   layer temperature (ta)                            k       m*ndim*np
!   layer specific humidity (wa)                      g/g     m*ndim*np
!   layer ozone mixing ratio by mass (oa)             g/g     m*ndim*np
!   surface temperature (ts)                          k       m*ndim
!   co2 mixing ratio by volumn (co2)                  pppv     1
!   high                                                       1
!
! pre-computed tables used in table look-up for transmittance calculations:
!
!   c1 , c2, c3: for co2 (band 3)
!   o1 , o2, o3: for  o3 (band 5)
!   h11,h12,h13: for h2o (band 1)
!   h21,h22,h23: for h2o (band 2)
!   h71,h72,h73: for h2o (band 7)
!
!---- output parameters
!
!   net downward flux, all-sky   (flx)             w/m**2     m*ndim*(np+1)
!   net downward flux, clear-sky (flc)             w/m**2     m*ndim*(np+1)
!   sensitivity of net downward flux
!    to surface temperature (dfdts)             w/m**2/k   m*ndim*(np+1)
!   emission by the surface (st4)                  w/m**2     m*ndim
!
! Notes:
!
!   (1)  Water vapor continuum absorption is included in 540-1380 /cm.
!   (2)  Scattering by clouds is not included.
!   (3)  Clouds are assumed "gray" bodies.
!   (4)  The diffuse cloud transmission is computed to be exp(-1.66*taucl).
!   (5)  If there are no clouds, flx=flc.
!   (6)  plevel(1) is the pressure at the top of the model atmosphere, and
!     plevel(np+1) is the surface pressure.
!
!    ARPS note: pl was replaced by pa at scalar points (layers)
!
!   (7)  Downward flux is positive, and upward flux is negative.
!   (8)  dfdts is always negative because upward flux is defined as negative.
!   (9)  For questions and coding errors, please contact with Ming-Dah Chou,
!     Code 913, NASA/Goddard Space Flight Center, Greenbelt, MD 20771.
!     Phone: 301-286-4012, Fax: 301-286-1759,
!     e-mail: chou@climate.gsfc.nasa.gov
!
!-----parameters defining the size of the pre-computed tables for transmittance
!  calculations using table look-up.
!
!  "nx" is the number of intervals in pressure
!  "no" is the number of intervals in o3 amount
!  "nc" is the number of intervals in co2 amount
!  "nh" is the number of intervals in h2o amount
!  "nt" is the number of copies to be made from the pre-computed
!       transmittance tables to reduce "memory-bank conflict"
!       in parallel machines and, hence, enhancing the speed of
!       computations using table look-up.
!       If such advantage does not exist, "nt" can be set to 1.
!***************************************************************************
!
!fpp$ expand (expmn)
!dir$ inline always expmn
!*$*  inline routine (expmn)
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: idim,jdim

  INTEGER :: nx,no,nc,nh,nt,m,n,np
  INTEGER :: i,j,k,ip,iw,it,ib,ik,iq,isb,k1,k2
  PARAMETER (nx=26,no=21,nc=24,nh=31,nt=7)

!---- input parameters ------

  REAL :: taucl(idim,jdim,np)
  REAL :: ccld(idim,jdim,np)
  REAL :: pl(idim,jdim,np+1)
  REAL :: ta(idim,jdim,np)
  REAL :: wa(idim,jdim,np)
  REAL :: oa(idim,jdim,np)
  REAL :: ts(idim,jdim)
  REAL :: co2

  LOGICAL :: high

!---- output parameters ------

  REAL :: flx(idim,jdim,np+1)
  REAL :: flc(idim,jdim,np+1)
  REAL :: dfdts(idim,jdim,np+1)
  REAL :: st4(idim,jdim)
!
!-----------------------------------------------------------------------
!
!  Temporary arrays
!
!-----------------------------------------------------------------------
!
  REAL :: fclr(m,n)
  REAL :: dbs(m,n)
  REAL :: trant(m,n)
  REAL :: th2o(m,n,6)
  REAL :: tcon(m,n,3)
  REAL :: tco2(m,n,6,2)

  REAL :: pa(m,n,np)
  REAL :: dt(m,n,np)
  REAL :: sh2o(m,n,np+1)
  REAL :: swpre(m,n,np+1)
  REAL :: swtem(m,n,np+1)
  REAL :: sco3(m,n,np+1)
  REAL :: scopre(m,n,np+1)
  REAL :: scotem(m,n,np+1)
  REAL :: dh2o(m,n,np)
  REAL :: dcont(m,n,np)
  REAL :: dco2(m,n,np)
  REAL :: do3(m,n,np)
  REAL :: flxu(m,n,np+1)
  REAL :: flxd(m,n,np+1)
  REAL :: clr(m,n,0:np+1)
  REAL :: blayer(m,n,0:np+1)

  REAL :: h2oexp(m,n,np,6)
  REAL :: conexp(m,n,np,3)

  REAL :: co2exp(m,n,np,6,2)
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  LOGICAL :: oznbnd
  LOGICAL :: co2bnd
  LOGICAL :: h2otbl
  LOGICAL :: conbnd

  REAL :: c1 (nx,nc,nt),c2 (nx,nc,nt),c3 (nx,nc,nt)
  REAL :: o1 (nx,no,nt),o2 (nx,no,nt),o3 (nx,no,nt)
  REAL :: h11(nx,nh,nt),h12(nx,nh,nt),h13(nx,nh,nt)
  REAL :: h21(nx,nh,nt),h22(nx,nh,nt),h23(nx,nh,nt)
  REAL :: h71(nx,nh,nt),h72(nx,nh,nt),h73(nx,nh,nt)

  REAL :: xx, w1,p1,dwe,dpe

!---- static data -----

  REAL :: cb(5,8)

!-----the following coefficients (table 2 of chou and suarez, 1995)
!  are for computing spectrally integtrated planck fluxes of
!  the 8 bands using eq. (22)

  DATA cb/                                                              &
      -2.6844E-1,-8.8994E-2, 1.5676E-3,-2.9349E-6, 2.2233E-9,           &
      3.7315E+1,-7.4758E-1, 4.6151E-3,-6.3260E-6, 3.5647E-9,            &
      3.7187E+1,-3.9085E-1,-6.1072E-4, 1.4534E-5,-1.6863E-8,            &
      -4.1928E+1, 1.0027E+0,-8.5789E-3, 2.9199E-5,-2.5654E-8,           &
      -4.9163E+1, 9.8457E-1,-7.0968E-3, 2.0478E-5,-1.5514E-8,           &
      -1.0345E+2, 1.8636E+0,-1.1753E-2, 2.7864E-5,-1.1998E-8,           &
      -6.9233E+0,-1.5878E-1, 3.9160E-3,-2.4496E-5, 4.9301E-8,           &
      1.1483E+2,-2.2376E+0, 1.6394E-2,-5.3672E-5, 6.6456E-8/

!-----copy tables to enhance the speed of co2 (band 3), o3 (band5),
!  and h2o (bands 1, 2, and 7 only) transmission calculations
!  using table look-up.

  LOGICAL :: first
  SAVE first
  DATA first /.true./
!
!-----------------------------------------------------------------------
!
!  Functions:
!
!-----------------------------------------------------------------------
!
  REAL :: expmn
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
!  include "h2o.tran3"
!  include "co2.tran3"
!  include "o3.tran3"

  COMMON /radtab001/ h11,h12,h13,h21,h22,h23,h71,h72,h73
  COMMON /radtab002/ c1,c2,c3
  COMMON /radtab003/ o1,o2,o3
!
!-----------------------------------------------------------------------
!
!  Save variables:
!
!-----------------------------------------------------------------------
!
!  save c1,c2,c3,o1,o2,o3
!  save h11,h12,h13,h21,h22,h23,h71,h72,h73
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  IF (first) THEN

!-----tables co2 and h2o are only used with 'high' option

    IF (high) THEN

      DO iw=1,nh
        DO ip=1,nx
          h11(ip,iw,1)=1.0-h11(ip,iw,1)
          h21(ip,iw,1)=1.0-h21(ip,iw,1)
          h71(ip,iw,1)=1.0-h71(ip,iw,1)
        END DO
      END DO

      DO iw=1,nc
        DO ip=1,nx
          c1(ip,iw,1)=1.0-c1(ip,iw,1)
        END DO
      END DO

!-----tables are replicated to avoid memory bank conflicts

      DO it=2,nt
        DO iw=1,nc
          DO ip=1,nx
            c1 (ip,iw,it)= c1(ip,iw,1)
            c2 (ip,iw,it)= c2(ip,iw,1)
            c3 (ip,iw,it)= c3(ip,iw,1)
          END DO
        END DO
        DO iw=1,nh
          DO ip=1,nx
            h11(ip,iw,it)=h11(ip,iw,1)
            h12(ip,iw,it)=h12(ip,iw,1)
            h13(ip,iw,it)=h13(ip,iw,1)
            h21(ip,iw,it)=h21(ip,iw,1)
            h22(ip,iw,it)=h22(ip,iw,1)
            h23(ip,iw,it)=h23(ip,iw,1)
            h71(ip,iw,it)=h71(ip,iw,1)
            h72(ip,iw,it)=h72(ip,iw,1)
            h73(ip,iw,it)=h73(ip,iw,1)
          END DO
        END DO
      END DO

    END IF

!-----always use table look-up for ozone transmittance

    DO iw=1,no
      DO ip=1,nx
        o1(ip,iw,1)=1.0-o1(ip,iw,1)
      END DO
    END DO

    DO it=2,nt
      DO iw=1,no
        DO ip=1,nx
          o1 (ip,iw,it)= o1(ip,iw,1)
          o2 (ip,iw,it)= o2(ip,iw,1)
          o3 (ip,iw,it)= o3(ip,iw,1)
        END DO
      END DO
    END DO

    first=.false.

  END IF

!-----compute layer pressure (pa) and layer temperature minus 250K (dt)

  DO k=1,np
    DO j=1,n
      DO i=1,m
        dt(i,j,k)=ta(i,j,k)-250.0
        pa(i,j,k)=0.5*(pl(i,j,k)+pl(i,j,k+1))
      END DO
    END DO
  END DO

!-----compute layer absorber amount

!  dh2o : water vapor amount (g/cm**2)
!  dcont: scaled water vapor amount for continuum absorption (g/cm**2)
!  dco2 : co2 amount (cm-atm)stp
!  do3  : o3 amount (cm-atm)stp
!  the factor 1.02 is equal to 1000/980
!  factors 789 and 476 are for unit conversion
!  the factor 0.001618 is equal to 1.02/(.622*1013.25)
!  the factor 6.081 is equal to 1800/296

  DO k=1,np
    DO j=1,n
      DO i=1,m
        dh2o(i,j,k) = 1.02*wa(i,j,k)*(pl(i,j,k+1)-pl(i,j,k))+1.e-10
        dco2(i,j,k) = 789.*co2*(pl(i,j,k+1)-pl(i,j,k))+1.e-10
        do3 (i,j,k) = 476.0*oa(i,j,k)*(pl(i,j,k+1)-pl(i,j,k))+1.e-10

!-----compute scaled water vapor amount for h2o continuum absorption
!  following eq. (43).

        xx=pa(i,j,k)*0.001618*wa(i,j,k)*wa(i,j,k)                       &
            *(pl(i,j,k+1)-pl(i,j,k))
        dcont(i,j,k) = xx*expmn(1800./ta(i,j,k)-6.081)+1.e-10

!-----compute effective cloud-free fraction, clr, for each layer.
!  the cloud diffuse transmittance is approximated by using a
!  diffusivity factor of 1.66.

        clr(i,j,k)=1.0-(ccld(i,j,k)*(1.-expmn(-1.66*taucl(i,j,k))))

      END DO
    END DO
  END DO

!-----compute column-integrated h2o amoumt, h2o-weighted pressure
!  and temperature.  it follows eqs. (37) and (38).

  IF (high) THEN

    CALL column(m,n,np,pa,dt,dh2o,sh2o,swpre,swtem)

  END IF

!-----the surface (with an index np+1) is treated as a layer filled with
!  black clouds.

  DO j=1,n
    DO i=1,m
      clr(i,j,0)    = 1.0
      clr(i,j,np+1) = 0.0
      st4(i,j)      = 0.0
    END DO
  END DO

!-----initialize fluxes

  DO k=1,np+1
    DO j=1,n
      DO i=1,m
        flx(i,j,k)  = 0.0
        flc(i,j,k)  = 0.0
        dfdts(i,j,k)= 0.0
        flxu(i,j,k) = 0.0
        flxd(i,j,k) = 0.0
      END DO
    END DO
  END DO

!-----integration over spectral bands

  DO ib=1,8

!-----if h2otbl, compute h2o (line) transmittance using table look-up.
!  if conbnd, compute h2o (continuum) transmittance in bands 3, 4, 5 and 6.
!  if co2bnd, compute co2 transmittance in band 3.
!  if oznbnd, compute  o3 transmittance in band 5.

    h2otbl=high.AND.(ib == 1.OR.ib == 2.OR.ib == 7)
    conbnd=ib >= 3.AND.ib <= 6
    co2bnd=ib == 3
    oznbnd=ib == 5

!-----blayer is the spectrally integrated planck flux of the mean layer
!  temperature derived from eq. (22)
!  the fitting for the planck flux is valid in the range 160-345 K.

    DO k=1,np
      DO j=1,n
        DO i=1,m
          blayer(i,j,k)=ta(i,j,k)*(ta(i,j,k)*(ta(i,j,k)                 &
                       *(ta(i,j,k)*cb(5,ib)+cb(4,ib))+cb(3,ib))         &
                       +cb(2,ib))+cb(1,ib)
        END DO
      END DO
    END DO

!-----the earth's surface, with an index "np+1", is treated as a layer

    DO j=1,n
      DO i=1,m
        blayer(i,j,0)   = 0.0
        blayer(i,j,np+1)=ts(i,j)*(ts(i,j)*(ts(i,j)                      &
                        *(ts(i,j)*cb(5,ib)+cb(4,ib))+cb(3,ib))          &
                        +cb(2,ib))+cb(1,ib)

!-----dbs is the derivative of the surface planck flux with respect to
!  surface temperature (eq. 59).

        dbs(i,j)=ts(i,j)*(ts(i,j)*(ts(i,j)                              &
                *4.*cb(5,ib)+3.*cb(4,ib))+2.*cb(3,ib))+cb(2,ib)

      END DO
    END DO

!-----compute column-integrated absorber amoumt, absorber-weighted
!  pressure and temperature for co2 (band 3) and o3 (band 5).
!  it follows eqs. (37) and (38).

!-----this is in the band loop to save storage

    IF( high .AND. co2bnd) THEN

      CALL column(m,n,np,pa,dt,dco2,sco3,scopre,scotem)

    END IF

    IF(oznbnd) THEN

      CALL column(m,n,np,pa,dt,do3,sco3,scopre,scotem)

    END IF

!-----compute the exponential terms (eq. 32) at each layer for
!  water vapor line absorption when k-distribution is used

    IF( .NOT. h2otbl) THEN

      CALL h2oexps(ib,m,n,np,dh2o,pa,dt,h2oexp)

    END IF

!-----compute the exponential terms (eq. 46) at each layer for
!  water vapor continuum absorption

    IF( conbnd) THEN

      CALL conexps(ib,m,n,np,dcont,conexp)

    END IF


!-----compute the  exponential terms (eq. 32) at each layer for
!  co2 absorption

    IF( .NOT.high .AND. co2bnd) THEN

      CALL co2exps(m,n,np,dco2,pa,dt,co2exp)

    END IF

!-----compute transmittances for regions between levels k1 and k2
!  and update the fluxes at the two levels.

    DO k1=1,np

!-----initialize fclr, th2o, tcon, and tco2

      DO j=1,n
        DO i=1,m
          fclr(i,j)=1.0
        END DO
      END DO

!-----for h2o line absorption

      IF(.NOT. h2otbl) THEN
        DO ik=1,6
          DO j=1,n
            DO i=1,m
              th2o(i,j,ik)=1.0
            END DO
          END DO
        END DO
      END IF

!-----for h2o continuum absorption

      IF (conbnd) THEN
        DO iq=1,3
          DO j=1,n
            DO i=1,m
              tcon(i,j,iq)=1.0
            END DO
          END DO
        END DO
      END IF

!-----for co2 absorption when using k-distribution method.
!  band 3 is divided into 3 sub-bands, but sub-bands 3a and 3c
!  are combined in computing the co2 transmittance.

      IF (.NOT. high .AND. co2bnd) THEN
        DO isb=1,2
          DO ik=1,6
            DO j=1,n
              DO i=1,m
                tco2(i,j,ik,isb)=1.0
              END DO
            END DO
          END DO
        END DO
      END IF

!-----loop over the bottom level of the region (k2)

      DO k2=k1+1,np+1

        DO j=1,n
          DO i=1,m
            trant(i,j)=1.0
          END DO
        END DO

        IF(h2otbl) THEN

          w1=-8.0
          p1=-2.0
          dwe=0.3
          dpe=0.2

!-----compute water vapor transmittance using table look-up

          IF (ib == 1 ) THEN

            CALL tablup(k1,k2,m,n,np,nx,nh,nt,sh2o,swpre,swtem,         &
                        w1,p1,dwe,dpe,h11,h12,h13,trant)

          END IF
          IF (ib == 2 ) THEN

            CALL tablup(k1,k2,m,n,np,nx,nh,nt,sh2o,swpre,swtem,         &
                        w1,p1,dwe,dpe,h21,h22,h23,trant)

          END IF
          IF (ib == 7 ) THEN

            CALL tablup(k1,k2,m,n,np,nx,nh,nt,sh2o,swpre,swtem,         &
                        w1,p1,dwe,dpe,h71,h72,h73,trant)

          END IF

        ELSE

!-----compute water vapor transmittance using k-distribution.

          CALL wvkdis(ib,m,n,np,k2-1,h2oexp,conexp,th2o,tcon,trant)

        END IF

        IF(co2bnd) THEN

          IF( high ) THEN

!-----compute co2 transmittance using table look-up method

            w1=-4.0
            p1=-2.0
            dwe=0.3
            dpe=0.2

            CALL tablup(k1,k2,m,n,np,nx,nc,nt,sco3,scopre,scotem,       &
                        w1,p1,dwe,dpe,c1,c2,c3,trant)

          ELSE

!-----compute co2 transmittance using k-distribution method

            CALL co2kdis(m,n,np,k2-1,co2exp,tco2,trant)

          END IF

        END IF

!-----compute o3 transmittance using table look-up

        IF (oznbnd) THEN

          w1=-6.0
          p1=-2.0
          dwe=0.3
          dpe=0.2

          CALL tablup(k1,k2,m,n,np,nx,no,nt,sco3,scopre,scotem,         &
                      w1,p1,dwe,dpe,o1,o2,o3,trant)

        END IF

!-----fclr is the clear line-of-sight between levels k1 and k2.
!  in computing fclr, clouds are assumed randomly overlapped
!  using eq. (10).

        DO j=1,n
          DO i=1,m
            fclr(i,j) = fclr(i,j)*clr(i,j,k2-1)
          END DO
        END DO

!-----compute upward and downward fluxes


!-----add "boundary" terms to the net downward flux.
!  these are the first terms on the right-hand-side of
!  eqs. (56a) and (56b).
!  downward fluxes are positive.

        IF (k2 == k1+1) THEN
          DO j=1,n
            DO i=1,m
              flc(i,j,k1)=flc(i,j,k1)-blayer(i,j,k1)
              flc(i,j,k2)=flc(i,j,k2)+blayer(i,j,k1)
            END DO
          END DO
        END IF

!-----add flux components involving the four layers above and below
!  the levels k1 and k2.  it follows eqs. (56a) and (56b).

        DO j=1,n
          DO i=1,m
            xx=trant(i,j)*(blayer(i,j,k2-1)-blayer(i,j,k2))
            flc(i,j,k1) =flc(i,j,k1)+xx
            xx=trant(i,j)*(blayer(i,j,k1-1)-blayer(i,j,k1))
            flc(i,j,k2) =flc(i,j,k2)+xx
          END DO
        END DO

!-----compute upward and downward fluxes for all-sky situation

        IF (k2 == k1+1) THEN
          DO j=1,n
            DO i=1,m
              flxu(i,j,k1)=flxu(i,j,k1)-blayer(i,j,k1)
              flxd(i,j,k2)=flxd(i,j,k2)+blayer(i,j,k1)
            END DO
          END DO
        END IF

        DO j=1,n
          DO i=1,m
            xx=trant(i,j)*(blayer(i,j,k2-1)-blayer(i,j,k2))
            flxu(i,j,k1) =flxu(i,j,k1)+xx*fclr(i,j)
            xx=trant(i,j)*(blayer(i,j,k1-1)-blayer(i,j,k1))
            flxd(i,j,k2) =flxd(i,j,k2)+xx*fclr(i,j)
          END DO
        END DO


      END DO

!-----compute the partial derivative of fluxes with respect to
!  surface temperature (eq. 59).

      DO j=1,n
        DO i=1,m
          dfdts(i,j,k1) =dfdts(i,j,k1)-dbs(i,j)*trant(i,j)*fclr(i,j)
        END DO
      END DO

    END DO

!-----add contribution from the surface to the flux terms at the surface.

    DO j=1,n
      DO i=1,m
        dfdts(i,j,np+1) =dfdts(i,j,np+1)-dbs(i,j)
      END DO
    END DO

    DO j=1,n
      DO i=1,m
        flc(i,j,np+1)=flc(i,j,np+1)-blayer(i,j,np+1)
        flxu(i,j,np+1)=flxu(i,j,np+1)-blayer(i,j,np+1)
        st4(i,j)=st4(i,j)-blayer(i,j,np+1)
      END DO
    END DO


!  write(7,3211) ib, flxd(1,1,52),flxu(1,1,52)
!  write(7,3211) ib, flxd(1,1,np+1),flxu(1,1,np+1)
!    3211 FORMAT ('ib, fluxd, fluxu=', i3,2F12.3)

  END DO

  DO k=1,np+1
    DO j=1,n
      DO i=1,m
        flx(i,j,k)=flxd(i,j,k)+flxu(i,j,k)
      END DO
    END DO
  END DO

  RETURN
END SUBROUTINE irrad
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE H2OEXPS                    ######
!######                                                      ######
!######                     Developed by                     ######
!######                                                      ######
!######    Goddard Cumulus Ensemble Modeling Group, NASA     ######
!######                                                      ######
!######     Center for Analysis and Prediction of Storms     ######
!######               University of Oklahoma                 ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE h2oexps(ib,m,n,np,dh2o,pa,dt,h2oexp) 1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate exponentials for water vapor line absorption in
!  individual layers.
!
!-----------------------------------------------------------------------
!
!
!  AUTHOR: (a) Radiative Transfer Model: M.-D. Chou and M. Suarez
!          (b) Cloud Optics:Tao, Lang, Simpson, Sui, Ferrier and
!              Chou (1996)
!
!  MODIFICATION:
!
!  03/11/1996 (Yuhe Liu)
!  Formatted code to ARPS standard format
!
!-----------------------------------------------------------------------
!
!  ORIGINAL COMMENTS:
!
!   compute exponentials for water vapor line absorption
!   in individual layers.
!
!---- input parameters
!  spectral band (ib)
!  number of grid intervals in zonal direction (m)
!  number of grid intervals in meridional direction (n)
!  number of layers (np)
!  layer water vapor amount for line absorption (dh2o)
!  layer pressure (pa)
!  layer temperature minus 250K (dt)
!
!---- output parameters
!  6 exponentials for each layer  (h2oexp)
!
!**********************************************************************
!
!fpp$ expand (expmn)
!dir$ inline always expmn
!*$*  inline routine (expmn)
!
!-----------------------------------------------------------------------
!
!  Variable declarations
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: ib,m,n,np

!---- input parameters ------

  REAL :: dh2o(m,n,np)
  REAL :: pa(m,n,np)
  REAL :: dt(m,n,np)

!---- output parameters -----

  REAL :: h2oexp(m,n,np,6)

!---- static data -----

  INTEGER :: mw(8)

  REAL :: xkw(8)
  REAL :: aw(8)
  REAL :: bw(8)

!---- temporary arrays -----

  REAL :: xh

!---- local misc. variables

  INTEGER :: i,j,k,ik

!-----xkw  are the absorption coefficients for the first
!  k-distribution function due to water vapor line absorption
!  (tables 4 and 7).  units are cm**2/g

  DATA xkw / 29.55  , 4.167E-1, 1.328E-2, 5.250E-4,                     &
              5.25E-4, 2.340E-3, 1.320E-0, 5.250E-4/

!-----mw are the ratios between neighboring absorption coefficients
!  for water vapor line absorption (tables 4 and 7).

  DATA mw /6,6,8,6,6,8,6,16/

!-----aw and bw (table 3) are the coefficients for temperature scaling
!  in eq. (25).

  DATA aw/ 0.0021, 0.0140, 0.0167, 0.0302,                              &
           0.0307, 0.0154, 0.0008, 0.0096/
  DATA bw/ -1.01E-5, 5.57E-5, 8.54E-5, 2.96E-4,                         &
            2.86E-4, 7.53E-5,-3.52E-6, 1.64E-5/

!-----expmn is an external function

  REAL :: expmn
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

!**********************************************************************
!    note that the 3 sub-bands in band 3 use the same set of xkw, aw,
!    and bw.  therefore, h2oexp for these sub-bands are identical.
!**********************************************************************

  DO k=1,np
    DO j=1,n
      DO i=1,m

!-----xh is   the scaled water vapor amount for line absorption
!  computed from (27).

        xh = dh2o(i,j,k)*(pa(i,j,k)*0.002)                              &
            * ( 1.+(aw(ib)+bw(ib)* dt(i,j,k))*dt(i,j,k) )

!-----h2oexp is the water vapor transmittance of the layer (k2-1)
!  due to line absorption

        h2oexp(i,j,k,1) = expmn(-xh*xkw(ib))

      END DO
    END DO
  END DO

  DO ik=2,6

    IF(mw(ib) == 6) THEN

      DO k=1,np
        DO j=1,n
          DO i=1,m
            xh = h2oexp(i,j,k,ik-1)*h2oexp(i,j,k,ik-1)
            h2oexp(i,j,k,ik) = xh*xh*xh
          END DO
        END DO
      END DO

    ELSE IF(mw(ib) == 8) THEN

      DO k=1,np
        DO j=1,n
          DO i=1,m
            xh = h2oexp(i,j,k,ik-1)*h2oexp(i,j,k,ik-1)
            xh = xh*xh
            h2oexp(i,j,k,ik) = xh*xh
          END DO
        END DO
      END DO

    ELSE

      DO k=1,np
        DO j=1,n
          DO i=1,m
            xh = h2oexp(i,j,k,ik-1)*h2oexp(i,j,k,ik-1)
            xh = xh*xh
            xh = xh*xh
            h2oexp(i,j,k,ik) = xh*xh
          END DO
        END DO
      END DO

    END IF
  END DO

  RETURN
END SUBROUTINE h2oexps
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE CONEXPS                    ######
!######                                                      ######
!######                     Developed by                     ######
!######                                                      ######
!######    Goddard Cumulus Ensemble Modeling Group, NASA     ######
!######                                                      ######
!######     Center for Analysis and Prediction of Storms     ######
!######               University of Oklahoma                 ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE conexps(ib,m,n,np,dcont,conexp) 1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate exponentials for continuum absorption in individual
!  layers.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: (a) Radiative Transfer Model: M.-D. Chou and M. Suarez
!          (b) Cloud Optics:Tao, Lang, Simpson, Sui, Ferrier and
!              Chou (1996)
!
!  MODIFICATION HISTORY:
!
!  03/15/1996 (Yuhe Liu)
!  Adopted the original code and formatted it in accordance with the
!  ARPS coding standard.
!
!-----------------------------------------------------------------------
!
!  ORIGINAL COMMENTS:
!
!**********************************************************************
!   compute exponentials for continuum absorption in individual layers.
!
!---- input parameters
!  spectral band (ib)
!  number of grid intervals in zonal direction (m)
!  number of grid intervals in meridional direction (n)
!  number of layers (np)
!  layer scaled water vapor amount for continuum absorption (dcont)
!
!---- output parameters
!  1 or 3 exponentials for each layer (conexp)
!
!**********************************************************************
!
!fpp$ expand (expmn)
!dir$ inline always expmn
!*$*  inline routine (expmn)
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: ib,m,n,np,i,j,k,iq

!---- input parameters ------

  REAL :: dcont(m,n,np)

!---- updated parameters -----

  REAL :: conexp(m,n,np,3)

!---- static data -----

  INTEGER :: NE(8)
  REAL :: xke(8)


!-----xke are the absorption coefficients for the first
!  k-distribution function due to water vapor continuum absorption
!  (table 6).  units are cm**2/g

  DATA xke /  0.00,   0.00,   27.40,   15.8,                            &
              9.40,   7.75,     0.0,    0.0/


!-----ne is the number of terms in computing water vapor
!  continuum transmittance (Table 6).
!  band 3 is divided into 3 sub-bands.

  DATA NE /0,0,3,1,1,1,0,0/
!
!-----------------------------------------------------------------------
!
!  Functions:
!
!-----------------------------------------------------------------------
!

!-----expmn is an external function
  REAL :: expmn
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO k=1,np
    DO j=1,n
      DO i=1,m
        conexp(i,j,k,1) = expmn(-dcont(i,j,k)*xke(ib))
      END DO
    END DO
  END DO

  IF (ib == 3) THEN

!-----the absorption coefficients for sub-bands 3b (iq=2) and 3a (iq=3)
!  are, respectively, double and quadruple that for sub-band 3c (iq=1)
!  (table 6).

    DO iq=2,3
      DO k=1,np
        DO j=1,n
          DO i=1,m
            conexp(i,j,k,iq) = conexp(i,j,k,iq-1) *conexp(i,j,k,iq-1)
          END DO
        END DO
      END DO
    END DO

  END IF

  RETURN
END SUBROUTINE conexps
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE CO2EXPS                    ######
!######                                                      ######
!######                     Developed by                     ######
!######                                                      ######
!######    Goddard Cumulus Ensemble Modeling Group, NASA     ######
!######                                                      ######
!######     Center for Analysis and Prediction of Storms     ######
!######               University of Oklahoma                 ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE co2exps(m,n,np,dco2,pa,dt,co2exp) 1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate co2 exponentials for individual layers.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: (a) Radiative Transfer Model: M.-D. Chou and M. Suarez
!          (b) Cloud Optics:Tao, Lang, Simpson, Sui, Ferrier and
!              Chou (1996)
!
!  MODIFICATION HISTORY:
!
!  03/15/1996 (Yuhe Liu)
!  Adopted the original code and formatted it in accordance with the
!  ARPS coding standard.
!
!-----------------------------------------------------------------------
!
!  ORIGINAL COMMENTS:
!
!
!**********************************************************************
!   compute co2 exponentials for individual layers.
!
!---- input parameters
!  number of grid intervals in zonal direction (m)
!  number of grid intervals in meridional direction (n)
!  number of layers (np)
!  layer co2 amount (dco2)
!  layer pressure (pa)
!  layer temperature minus 250K (dt)
!
!---- output parameters
!  6 exponentials for each layer (co2exp)
!**********************************************************************
!
!fpp$ expand (expmn)
!dir$ inline always expmn
!*$*  inline routine (expmn)
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: m,n,np,i,j,k

!---- input parameters -----

  REAL :: dco2(m,n,np)
  REAL :: pa(m,n,np)
  REAL :: dt(m,n,np)

!---- output parameters -----

  REAL :: co2exp(m,n,np,6,2)

!---- static data -----

  REAL :: xkc(2),ac(2),bc(2),pm(2),prc(2)

!---- temporary arrays -----

  REAL :: xc

!-----xkc is the absorption coefficients for the
!  first k-distribution function due to co2 (table 7).
!  units are 1/(cm-atm)stp.

  DATA xkc/2.656E-5,2.656E-3/

!-----parameters (table 3) for computing the scaled co2 amount
!  using (27).

  DATA prc/  300.0,   30.0/
  DATA pm /    0.5,   0.85/
  DATA ac / 0.0182, 0.0042/
  DATA bc /1.07E-4,2.00E-5/
!
!-----------------------------------------------------------------------
!
!  Functions:
!
!-----------------------------------------------------------------------
!

!-----function expmn

  REAL :: expmn
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO k=1,np
    DO j=1,n
      DO i=1,m

!-----compute the scaled co2 amount from eq. (27) for band-wings
!  (sub-bands 3a and 3c).

        xc = dco2(i,j,k)*(pa(i,j,k)/prc(1))**pm(1)                      &
                *(1.+(ac(1)+bc(1)*dt(i,j,k))*dt(i,j,k))

!-----six exponential by powers of 8 (table 7).

        co2exp(i,j,k,1,1)=expmn(-xc*xkc(1))

        xc=co2exp(i,j,k,1,1)*co2exp(i,j,k,1,1)
        xc=xc*xc
        co2exp(i,j,k,2,1)=xc*xc

        xc=co2exp(i,j,k,2,1)*co2exp(i,j,k,2,1)
        xc=xc*xc
        co2exp(i,j,k,3,1)=xc*xc

        xc=co2exp(i,j,k,3,1)*co2exp(i,j,k,3,1)
        xc=xc*xc
        co2exp(i,j,k,4,1)=xc*xc

        xc=co2exp(i,j,k,4,1)*co2exp(i,j,k,4,1)
        xc=xc*xc
        co2exp(i,j,k,5,1)=xc*xc

        xc=co2exp(i,j,k,5,1)*co2exp(i,j,k,5,1)
        xc=xc*xc
        co2exp(i,j,k,6,1)=xc*xc

!-----compute the scaled co2 amount from eq. (27) for band-center
!  region (sub-band 3b).

        xc = dco2(i,j,k)*(pa(i,j,k)/prc(2))**pm(2)                      &
                *(1.+(ac(2)+bc(2)*dt(i,j,k))*dt(i,j,k))

        co2exp(i,j,k,1,2)=expmn(-xc*xkc(2))

        xc=co2exp(i,j,k,1,2)*co2exp(i,j,k,1,2)
        xc=xc*xc
        co2exp(i,j,k,2,2)=xc*xc

        xc=co2exp(i,j,k,2,2)*co2exp(i,j,k,2,2)
        xc=xc*xc
        co2exp(i,j,k,3,2)=xc*xc

        xc=co2exp(i,j,k,3,2)*co2exp(i,j,k,3,2)
        xc=xc*xc
        co2exp(i,j,k,4,2)=xc*xc

        xc=co2exp(i,j,k,4,2)*co2exp(i,j,k,4,2)
        xc=xc*xc
        co2exp(i,j,k,5,2)=xc*xc

        xc=co2exp(i,j,k,5,2)*co2exp(i,j,k,5,2)
        xc=xc*xc
        co2exp(i,j,k,6,2)=xc*xc

      END DO
    END DO
  END DO

  RETURN
END SUBROUTINE co2exps
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE COLUMN                     ######
!######                                                      ######
!######                     Developed by                     ######
!######                                                      ######
!######    Goddard Cumulus Ensemble Modeling Group, NASA     ######
!######                                                      ######
!######     Center for Analysis and Prediction of Storms     ######
!######               University of Oklahoma                 ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE column(m,n,np,pa,dt,sabs0,sabs,spre,stem) 3
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate column-integrated (from top of the model atmosphere)
!  absorber amount, absorber-weighted pressure and temperature.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: (a) Radiative Transfer Model: M.-D. Chou and M. Suarez
!          (b) Cloud Optics:Tao, Lang, Simpson, Sui, Ferrier and
!              Chou (1996)
!
!  MODIFICATION HISTORY:
!
!  03/15/1996 (Yuhe Liu)
!  Adopted the original code and formatted it in accordance with the
!  ARPS coding standard.
!
!-----------------------------------------------------------------------
!
!  ORIGINAL COMMENTS:
!
!**************************************************************************
!-----compute column-integrated (from top of the model atmosphere)
!  absorber amount (sabs), absorber-weighted pressure (spre) and
!  temperature (stem).
!  computations of spre and stem follows eqs. (37) and (38).
!
!--- input parameters
!   number of soundings in zonal direction (m)
!   number of soundings in meridional direction (n)
!   number of atmospheric layers (np)
!   layer pressure (pa)
!   layer temperature minus 250K (dt)
!   layer absorber amount (sabs0)
!
!--- output parameters
!   column-integrated absorber amount (sabs)
!   column absorber-weighted pressure (spre)
!   column absorber-weighted temperature (stem)
!
!--- units of pa and dt are mb and k, respectively.
!    units of sabs are g/cm**2 for water vapor and (cm-atm)stp for co2 and o3
!**************************************************************************
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: m,n,np,i,j,k

!---- input parameters -----

  REAL :: pa(m,n,np)
  REAL :: dt(m,n,np)
  REAL :: sabs0(m,n,np)

!---- output parameters -----

  REAL :: sabs(m,n,np+1)
  REAL :: spre(m,n,np+1)
  REAL :: stem(m,n,np+1)

!*********************************************************************
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO j=1,n
    DO i=1,m
      sabs(i,j,1)=0.0
      spre(i,j,1)=0.0
      stem(i,j,1)=0.0
    END DO
  END DO

  DO k=1,np
    DO j=1,n
      DO i=1,m
        sabs(i,j,k+1)=sabs(i,j,k)+sabs0(i,j,k)
        spre(i,j,k+1)=spre(i,j,k)+pa(i,j,k)*sabs0(i,j,k)
        stem(i,j,k+1)=stem(i,j,k)+dt(i,j,k)*sabs0(i,j,k)
      END DO
    END DO
  END DO

  RETURN
END SUBROUTINE column
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE TABLUP                     ######
!######                                                      ######
!######                     Developed by                     ######
!######                                                      ######
!######    Goddard Cumulus Ensemble Modeling Group, NASA     ######
!######                                                      ######
!######     Center for Analysis and Prediction of Storms     ######
!######               University of Oklahoma                 ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE tablup(k1,k2,m,n,np,nx,nh,nt,sabs,spre,stem,w1,p1,           & 5
           dwe,dpe,coef1,coef2,coef3,tran)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate water vapor, co2, and o3 transmittances using table
!  look-up. rlwopt = 0, or high=.false.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: (a) Radiative Transfer Model: M.-D. Chou and M. Suarez
!          (b) Cloud Optics:Tao, Lang, Simpson, Sui, Ferrier and
!              Chou (1996)
!
!  MODIFICATION HISTORY:
!
!  03/15/1996 (Yuhe Liu)
!  Adopted the original code and formatted it in accordance with the
!  ARPS coding standard.
!
!-----------------------------------------------------------------------
!
!  ORIGINAL COMMENTS:
!
!**********************************************************************
!   compute water vapor, co2, and o3 transmittances between levels k1 and k2
!   using table look-up for m x n soundings.
!
!   Calculations follow Eq. (40) of Chou and Suarez (1995)
!
!---- input ---------------------
!  indices for pressure levels (k1 and k2)
!  number of grid intervals in zonal direction (m)
!  number of grid intervals in meridional direction (n)
!  number of atmospheric layers (np)
!  number of pressure intervals in the table (nx)
!  number of absorber amount intervals in the table (nh)
!  number of tables copied (nt)
!  column-integrated absorber amount (sabs)
!  column absorber amount-weighted pressure (spre)
!  column absorber amount-weighted temperature (stem)
!  first value of absorber amount (log10) in the table (w1)
!  first value of pressure (log10) in the table (p1)
!  size of the interval of absorber amount (log10) in the table (dwe)
!  size of the interval of pressure (log10) in the table (dpe)
!  pre-computed coefficients (coef1, coef2, and coef3)
!
!---- updated ---------------------
!  transmittance (tran)
!
!  Note:
!   (1) units of sabs are g/cm**2 for water vapor and (cm-atm)stp for co2 and o3.
!   (2) units of spre and stem are, respectively, mb and K.
!   (3) there are nt identical copies of the tables (coef1, coef2, and
!    coef3).  the prupose of using the multiple copies of tables is
!    to increase the speed in parallel (vectorized) computations.
!    if such advantage does not exist, nt can be set to 1.
!
!**********************************************************************
!
!-----------------------------------------------------------------------
!
!  Variable declarations
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: k1,k2,m,n,np,nx,nh,nt,i,j

!---- input parameters -----

  REAL :: w1,p1,dwe,dpe
  REAL :: sabs(m,n,np+1)
  REAL :: spre(m,n,np+1)
  REAL :: stem(m,n,np+1)

  REAL :: coef1(nx,nh,nt)
  REAL :: coef2(nx,nh,nt)
  REAL :: coef3(nx,nh,nt)

!---- update parameter -----

  REAL :: tran(m,n)

!---- temporary variables -----

  REAL :: x1,x2,x3,we,pe,fw,fp,pa,pb,pc,ax,ba,bb,t1,ca,cb,t2
  INTEGER :: iw,ip,nn

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO j=1,n
    DO i=1,m

      nn=MOD(i,nt)+1

      x1=sabs(i,j,k2)-sabs(i,j,k1)
      x2=(spre(i,j,k2)-spre(i,j,k1))/x1
      x3=(stem(i,j,k2)-stem(i,j,k1))/x1

      we=(LOG10(x1)-w1)/dwe
      pe=(LOG10(x2)-p1)/dpe

      we=MAX(we,w1-2.*dwe)
      pe=MAX(pe,p1)

      iw=INT(we+1.5)
      ip=INT(pe+1.5)

      iw=MIN(iw,nh-1)
      iw=MAX(iw, 2)

      ip=MIN(ip,nx-1)
      ip=MAX(ip, 1)

      fw=we-FLOAT(iw-1)
      fp=pe-FLOAT(ip-1)

!-----linear interpolation in pressure

      pa = coef1(ip,iw-1,nn)*(1.-fp)+coef1(ip+1,iw-1,nn)*fp
      pb = coef1(ip,iw,  nn)*(1.-fp)+coef1(ip+1,iw,  nn)*fp
      pc = coef1(ip,iw+1,nn)*(1.-fp)+coef1(ip+1,iw+1,nn)*fp

!-----quadratic interpolation in absorber amount for coef1

      ax = (-pa*(1.-fw)+pc*(1.+fw)) *fw*0.5 + pb*(1.-fw*fw)

!-----linear interpolation in absorber amount for coef2 and coef3

      ba = coef2(ip,iw,  nn)*(1.-fp)+coef2(ip+1,iw,  nn)*fp
      bb = coef2(ip,iw+1,nn)*(1.-fp)+coef2(ip+1,iw+1,nn)*fp
      t1 = ba*(1.-fw) + bb*fw

      ca = coef3(ip,iw,  nn)*(1.-fp)+coef3(ip+1,iw,  nn)*fp
      cb = coef3(ip,iw+1,nn)*(1.-fp)+coef3(ip+1,iw+1,nn)*fp
      t2 = ca*(1.-fw) + cb*fw

!-----update the total transmittance between levels k1 and k2

      tran(i,j)= (ax + (t1+t2*x3) * x3)*tran(i,j)

    END DO
  END DO

  RETURN
END SUBROUTINE tablup
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE WVKDIS                     ######
!######                                                      ######
!######                     Developed by                     ######
!######                                                      ######
!######    Goddard Cumulus Ensemble Modeling Group, NASA     ######
!######                                                      ######
!######     Center for Analysis and Prediction of Storms     ######
!######               University of Oklahoma                 ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE wvkdis(ib,m,n,np,k,h2oexp,conexp,th2o,tcon,tran) 1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate water vapor transmittance using the k-distribution
!  method.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: (a) Radiative Transfer Model: M.-D. Chou and M. Suarez
!          (b) Cloud Optics:Tao, Lang, Simpson, Sui, Ferrier and
!              Chou (1996)
!
!  MODIFICATION HISTORY:
!
!  03/15/1996 (Yuhe Liu)
!  Adopted the original code and formatted it in accordance with the
!  ARPS coding standard.
!
!-----------------------------------------------------------------------
!
!  ORIGINAL COMMENTS:
!
!**********************************************************************
!   compute water vapor transmittance between levels k1 and k2 for
!   m x n soundings using the k-distribution method.
!
!   computations follow eqs. (34), (46), (50) and (52).
!
!---- input parameters
!  spectral band (ib)
!  number of grid intervals in zonal direction (m)
!  number of grid intervals in meridional direction (n)
!  number of levels (np)
!  current level (k)
!  exponentials for line absorption (h2oexp)
!  exponentials for continuum absorption (conexp)
!
!---- updated parameters
!  transmittance between levels k1 and k2 due to
!    water vapor line absorption (th2o)
!  transmittance between levels k1 and k2 due to
!    water vapor continuum absorption (tcon)
!  total transmittance (tran)
!
!**********************************************************************
!
!-----------------------------------------------------------------------
!
!  Variable declarations
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: ib,m,n,np,k,i,j

!---- input parameters ------

  REAL :: conexp(m,n,np,3)
  REAL :: h2oexp(m,n,np,6)

!---- updated parameters -----

  REAL :: th2o(m,n,6)
  REAL :: tcon(m,n,3)
  REAL :: tran(m,n)

!---- static data -----

  INTEGER :: NE(8)
  REAL :: fkw(6,8),gkw(6,3)

!---- temporary variable -----

  REAL :: trnth2o

!-----fkw is the planck-weighted k-distribution function due to h2o
!  line absorption given in table 4 of Chou and Suarez (1995).
!  the k-distribution function for the third band, fkw(*,3), is not used

  DATA fkw / 0.2747,0.2717,0.2752,0.1177,0.0352,0.0255,                 &
             0.1521,0.3974,0.1778,0.1826,0.0374,0.0527,                 &
             6*1.00,                                                    &
             0.4654,0.2991,0.1343,0.0646,0.0226,0.0140,                 &
             0.5543,0.2723,0.1131,0.0443,0.0160,0.0000,                 &
             0.1846,0.2732,0.2353,0.1613,0.1146,0.0310,                 &
             0.0740,0.1636,0.4174,0.1783,0.1101,0.0566,                 &
             0.1437,0.2197,0.3185,0.2351,0.0647,0.0183/

!-----gkw is the planck-weighted k-distribution function due to h2o
!  line absorption in the 3 subbands (800-720,620-720,540-620 /cm)
!  of band 3 given in table 7.  Note that the order of the sub-bands
!  is reversed.

  DATA gkw/  0.1782,0.0593,0.0215,0.0068,0.0022,0.0000,                 &
             0.0923,0.1675,0.0923,0.0187,0.0178,0.0000,                 &
             0.0000,0.1083,0.1581,0.0455,0.0274,0.0041/

!-----ne is the number of terms used in each band to compute water vapor
!  continuum transmittance (table 6).

  DATA NE /0,0,3,1,1,1,0,0/
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

!-----tco2 are the six exp factors between levels k1 and k2
!  tran is the updated total transmittance between levels k1 and k2


!-----th2o is the 6 exp factors between levels k1 and k2 due to
!  h2o line absorption.

!-----tcon is the 3 exp factors between levels k1 and k2 due to
!  h2o continuum absorption.

!-----trnth2o is the total transmittance between levels k1 and k2 due
!  to both line and continuum absorption computed from eq. (52).

  DO j=1,n
    DO i=1,m
      th2o(i,j,1) = th2o(i,j,1)*h2oexp(i,j,k,1)
      th2o(i,j,2) = th2o(i,j,2)*h2oexp(i,j,k,2)
      th2o(i,j,3) = th2o(i,j,3)*h2oexp(i,j,k,3)
      th2o(i,j,4) = th2o(i,j,4)*h2oexp(i,j,k,4)
      th2o(i,j,5) = th2o(i,j,5)*h2oexp(i,j,k,5)
      th2o(i,j,6) = th2o(i,j,6)*h2oexp(i,j,k,6)
    END DO
  END DO


  IF (NE(ib) == 0) THEN


    DO j=1,n
      DO i=1,m

        trnth2o      =(fkw(1,ib)*th2o(i,j,1)                            &
                     + fkw(2,ib)*th2o(i,j,2)                            &
                     + fkw(3,ib)*th2o(i,j,3)                            &
                     + fkw(4,ib)*th2o(i,j,4)                            &
                     + fkw(5,ib)*th2o(i,j,5)                            &
                     + fkw(6,ib)*th2o(i,j,6))

        tran(i,j)=tran(i,j)*trnth2o

      END DO
    END DO

  ELSE IF (NE(ib) == 1) THEN


    DO j=1,n
      DO i=1,m

        tcon(i,j,1)= tcon(i,j,1)*conexp(i,j,k,1)

        trnth2o      =(fkw(1,ib)*th2o(i,j,1)                            &
                     + fkw(2,ib)*th2o(i,j,2)                            &
                     + fkw(3,ib)*th2o(i,j,3)                            &
                     + fkw(4,ib)*th2o(i,j,4)                            &
                     + fkw(5,ib)*th2o(i,j,5)                            &
                     + fkw(6,ib)*th2o(i,j,6))*tcon(i,j,1)

        tran(i,j)=tran(i,j)*trnth2o

      END DO
    END DO

  ELSE

    DO j=1,n
      DO i=1,m

        tcon(i,j,1)= tcon(i,j,1)*conexp(i,j,k,1)
        tcon(i,j,2)= tcon(i,j,2)*conexp(i,j,k,2)
        tcon(i,j,3)= tcon(i,j,3)*conexp(i,j,k,3)

        trnth2o      = (  gkw(1,1)*th2o(i,j,1)                          &
                        + gkw(2,1)*th2o(i,j,2)                          &
                        + gkw(3,1)*th2o(i,j,3)                          &
                        + gkw(4,1)*th2o(i,j,4)                          &
                        + gkw(5,1)*th2o(i,j,5)                          &
                        + gkw(6,1)*th2o(i,j,6) ) * tcon(i,j,1)          &
                     + (  gkw(1,2)*th2o(i,j,1)                          &
                        + gkw(2,2)*th2o(i,j,2)                          &
                        + gkw(3,2)*th2o(i,j,3)                          &
                        + gkw(4,2)*th2o(i,j,4)                          &
                        + gkw(5,2)*th2o(i,j,5)                          &
                        + gkw(6,2)*th2o(i,j,6) ) * tcon(i,j,2)          &
                     + (  gkw(1,3)*th2o(i,j,1)                          &
                        + gkw(2,3)*th2o(i,j,2)                          &
                        + gkw(3,3)*th2o(i,j,3)                          &
                        + gkw(4,3)*th2o(i,j,4)                          &
                        + gkw(5,3)*th2o(i,j,5)                          &
                        + gkw(6,3)*th2o(i,j,6) ) * tcon(i,j,3)

        tran(i,j)=tran(i,j)*trnth2o

      END DO
    END DO

  END IF

  RETURN
END SUBROUTINE wvkdis
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE CO2KDIS                    ######
!######                                                      ######
!######                     Developed by                     ######
!######                                                      ######
!######    Goddard Cumulus Ensemble Modeling Group, NASA     ######
!######                                                      ######
!######     Center for Analysis and Prediction of Storms     ######
!######               University of Oklahoma                 ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE co2kdis(m,n,np,k,co2exp,tco2,tran) 1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate co2 transmittances using the k-distribution method
!
!
!-----------------------------------------------------------------------
!
!  AUTHOR: (a) Radiative Transfer Model: M.-D. Chou and M. Suarez
!          (b) Cloud Optics:Tao, Lang, Simpson, Sui, Ferrier and
!              Chou (1996)
!
!  MODIFICATION HISTORY:
!
!  03/15/1996 (Yuhe Liu)
!  Adopted the original code and formatted it in accordance with the
!  ARPS coding standard.
!
!-----------------------------------------------------------------------
!
!  ORIGINAL COMMENTS:
!
!**********************************************************************
!   compute co2 transmittances between levels k1 and k2 for m x n soundings
!   using the k-distribution method with linear pressure scaling.
!
!   computations follow eq. (34).
!
!---- input parameters
!   number of grid intervals in zonal direction (m)
!   number of grid intervals in meridional direction (n)
!
!---- updated parameters
!   transmittance between levels k1 and k2 due to co2 absorption
!  for the various values of the absorption coefficient (tco2)
!   total transmittance (tran)
!
!**********************************************************************
!
!-----------------------------------------------------------------------
!
!  Variable declarations
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: m,n,np,k,i,j

!---- input parameters -----

  REAL :: co2exp(m,n,np,6,2)

!---- updated parameters -----

  REAL :: tco2(m,n,6,2)
  REAL :: tran(m,n)

!---- static data -----

  REAL :: gkc(6,2)

!---- temporary variable -----

  REAL :: xc

!-----gkc is the planck-weighted co2 k-distribution function
!  in the band-wing and band-center regions given in table 7.
!  for computing efficiency, sub-bands 3a and 3c are combined.

  DATA gkc/  0.1395,0.1407,0.1549,0.1357,0.0182,0.0220,                 &
             0.0766,0.1372,0.1189,0.0335,0.0169,0.0059/
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

!-----tco2 is the 6 exp factors between levels k1 and k2.
!  xc is the total co2 transmittance given by eq. (53).

  DO j=1,n
    DO i=1,m

!-----band-wings

      tco2(i,j,1,1)=tco2(i,j,1,1)*co2exp(i,j,k,1,1)
      xc=             gkc(1,1)*tco2(i,j,1,1)

      tco2(i,j,2,1)=tco2(i,j,2,1)*co2exp(i,j,k,2,1)
      xc=xc+gkc(2,1)*tco2(i,j,2,1)

      tco2(i,j,3,1)=tco2(i,j,3,1)*co2exp(i,j,k,3,1)
      xc=xc+gkc(3,1)*tco2(i,j,3,1)

      tco2(i,j,4,1)=tco2(i,j,4,1)*co2exp(i,j,k,4,1)
      xc=xc+gkc(4,1)*tco2(i,j,4,1)

      tco2(i,j,5,1)=tco2(i,j,5,1)*co2exp(i,j,k,5,1)
      xc=xc+gkc(5,1)*tco2(i,j,5,1)

      tco2(i,j,6,1)=tco2(i,j,6,1)*co2exp(i,j,k,6,1)
      xc=xc+gkc(6,1)*tco2(i,j,6,1)

!-----band-center region

      tco2(i,j,1,2)=tco2(i,j,1,2)*co2exp(i,j,k,1,2)
      xc=xc+gkc(1,2)*tco2(i,j,1,2)

      tco2(i,j,2,2)=tco2(i,j,2,2)*co2exp(i,j,k,2,2)
      xc=xc+gkc(2,2)*tco2(i,j,2,2)

      tco2(i,j,3,2)=tco2(i,j,3,2)*co2exp(i,j,k,3,2)
      xc=xc+gkc(3,2)*tco2(i,j,3,2)

      tco2(i,j,4,2)=tco2(i,j,4,2)*co2exp(i,j,k,4,2)
      xc=xc+gkc(4,2)*tco2(i,j,4,2)

      tco2(i,j,5,2)=tco2(i,j,5,2)*co2exp(i,j,k,5,2)
      xc=xc+gkc(5,2)*tco2(i,j,5,2)

      tco2(i,j,6,2)=tco2(i,j,6,2)*co2exp(i,j,k,6,2)
      xc=xc+gkc(6,2)*tco2(i,j,6,2)

      tran(i,j)=tran(i,j)*xc

    END DO
  END DO

  RETURN
END SUBROUTINE co2kdis