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


SUBROUTINE sorad(idim,jdim,m,n,np,                                      & 1,4
           pl,ta,wa,oa,co2,                                             &
           taucldi,taucldl,reffi,reffl,fcld,ict,icb,                    &
           taual,rsirbm,rsirdf,rsuvbm,rsuvdf,cosz,                      &
           flx,flc,fdirir,fdifir,fdirpar,fdifpar,                       &
           sdf,sclr,csm,cc,                                             &
           tauclb,tauclf,dp,wh,oh,scal,swh,so2,df,                      &
           tem2d1, tem2d2, tem2d3, tem2d4, tem2d5,                      &
           tem2d6, tem2d7, tem2d8, tem2d9, tem2d10,                     &
           tem2d11,tem2d12,tem2d13,tem2d14,tem2d15,                     &
           tem2d16,tem2d17,tem2d18,tem2d19,                             &
           tem3d1, tem3d2, tem3d3, tem3d4, tem3d5,                      &
           tem4d1, tem4d2, tem4d3, tem4d4, tem4d5,                      &
           tem5d1, tem5d2, tem5d3, tem5d4, tem5d5)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate solar fluxes due to the absoption by water vapor, ozone,
!  co2, o2, clouds, and aerosols and due to the scattering by clouds,
!  aerosols, and gases.
!
!-----------------------------------------------------------------------
!
!  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:
!
!****************  charney  /silo/z1mhy/new/solar.f  ***10/24/95******
!********************************************************************
!
! This routine computes solar fluxes due to the absoption by water
!  vapor, ozone, co2, o2, clouds, and aerosols and due to the
!  scattering by clouds, aerosols, and gases.
!
! This is a vectorized code. It computes the fluxes simultaneous for
!  (m x n) soundings, which is a subset of the (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.
!
! Ice and liquid cloud particles are allowed to co-exist in any of the
!  np layers. Two sets of cloud parameters are required as inputs, one
!  for ice paticles and the other for liquid particles.  These parameters
!  are optical thickness (taucld) and effective particle size (reff).
!
! If no information is available for reff, a default value of
!  10 micron for liquid water and 75 micron for ice can be used.
!
! Clouds are grouped into high, middle, and low clouds separated by the
!  level indices ict and icb.  For detail, see the subroutine cldscale.
!
!----- 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                   n/d        1
!        meridional direction (ndim)
!   number of atmospheric layers (np)                n/d        1
!   level pressure (pl)                              mb       m*ndim*(np+1)
!   layer temperature (ta)                           k        m*ndim*np
!   layer specific humidity (wa)                     gm/gm    m*ndim*np
!   layer ozone concentration (oa)                   gm/gm    m*ndim*np
!   co2 mixing ratio by volumn (co2)               parts/part   1
!   cloud optical thickness (taucld)                 n/d      m*ndim*np*2
!             index 1 for ice particles
!             index 2 for liquid drops
!   effective cloud-particle size (reff)           micrometer m*ndim*np*2
!             index 1 for ice particles
!             index 2 for liquid drops
!
!   ARPS notes: taucld were changed to taucldi and taucldl
!            reff   were changed to reffi   and reffl
!
!   cloud amount (fcld)                            fraction   m*ndim*np
!   level index separating high and middle           n/d        1
!             clouds (ict)
!   level index separating middle and low clouds     n/d        1
!             clouds (icb)
!   aerosol optical thickness (taual)                n/d      m*ndim*np
!   solar ir surface albedo for beam                fraction   m*ndim
!             radiation (rsirbm)
!   solar ir surface albedo for diffuse             fraction   m*ndim
!             radiation (rsirdf)
!   uv + par surface albedo for beam                     fraction   m*ndim
!             radiation (rsuvbm)
!   uv + par surface albedo for diffuse                  fraction   m*ndim
!             radiation (rsuvdf)
!   cosine of solar zenith angle (cosz)            n/d        m*ndim
!
!----- Output parameters
!
!   all-sky flux (downward minus upward) (flx)     fraction   m*ndim*(np+1)
!   clear-sky flux (downward minus upward) (flc)   fraction   m*ndim*(np+1)
!   all-sky direct downward ir (0.7-10 micron)
!             flux at the surface (fdirir)      fraction   m*ndim
!   all-sky diffuse downward ir flux at
!             the surface (fdifir)              fraction   m*ndim
!   all-sky direct downward par (0.4-0.7 micron)
!             flux at the surface (fdirpar)     fraction   m*ndim
!   all-sky diffuse downward par flux at
!             the surface (fdifpar)             fraction   m*ndim
!
!----- Notes:
!
!    (1) The unit of flux is fraction of the incoming solar radiation
!     at the top of the atmosphere.  Therefore, fluxes should
!     be equal to flux multiplied by the extra-terrestrial solar
!     flux and the cosine of solar zenith angle.
!    (2) Clouds and aerosols can be included in any layers by specifying
!     fcld(i,j,k), taucld(i,j,k,*) and taual(i,j,k), k=1,np.
!     For an atmosphere without clouds and aerosols,
!     set fcld(i,j,k)=taucld(i,j,k,*)=taual(i,j,k)=0.0.
!    (3) Aerosol single scattering albedos and asymmetry
!     factors are specified in the subroutines solir and soluv.
!    (4) pl(i,j,1) is the pressure at the top of the model, and
!     pl(i,j,np+1) is the surface pressure.
!
!    ARPS note: pl was replaced by pa at scalar points (layers)
!
!    (5) the pressure levels ict and icb correspond approximately
!     to 400 and 700 mb.
!
!**************************************************************************
!
!fpp$ expand (expmn)
!dir$ inline always expmn
!*$*  inline routine (expmn)
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: idim,jdim             ! ARPS nx, ny. NX is used in this
                                   ! subroutine for other purpose

  INTEGER :: m,n,np
  INTEGER :: ict,icb

  REAL :: pl(idim,jdim,np+1)       ! pressure at scalar points (layers)
  REAL :: ta(idim,jdim,np)
  REAL :: wa(idim,jdim,np)
  REAL :: oa(idim,jdim,np)
  REAL :: taucldi(idim,jdim,np)
  REAL :: taucldl(idim,jdim,np)
  REAL :: reffi(idim,jdim,np)
  REAL :: reffl(idim,jdim,np)
  REAL :: fcld(idim,jdim,np)
  REAL :: taual(idim,jdim,np)
  REAL :: rsirbm(idim,jdim)
  REAL :: rsirdf(idim,jdim)
  REAL :: rsuvbm(idim,jdim)
  REAL :: rsuvdf(idim,jdim)
  REAL :: cosz(idim,jdim)
  REAL :: co2

  REAL :: flx(idim,jdim,np+1)
  REAL :: flc(idim,jdim,np+1)
  REAL :: fdirir(idim,jdim)
  REAL :: fdifir(idim,jdim)
  REAL :: fdirpar(idim,jdim)
  REAL :: fdifpar(idim,jdim)
!
!-----------------------------------------------------------------------
!
!  Local temporary arrays
!
!-----------------------------------------------------------------------
!
  REAL :: sdf(m,n)
  REAL :: sclr(m,n)
  REAL :: csm(m,n)
  REAL :: cc(m,n,3)

  REAL :: tauclb(m,n,np)
  REAL :: tauclf(m,n,np)
  REAL :: dp(m,n,np)
  REAL :: wh(m,n,np)
  REAL :: oh(m,n,np)
  REAL :: scal(m,n,np)
  REAL :: swh(m,n,np+1)
  REAL :: so2(m,n,np+1)
  REAL :: df(m,n,np+1)

!-----temporary arrays used in subroutine cldflx

  REAL :: tem2d1 (m,n)
  REAL :: tem2d2 (m,n)
  REAL :: tem2d3 (m,n)
  REAL :: tem2d4 (m,n)
  REAL :: tem2d5 (m,n)
  REAL :: tem2d6 (m,n)
  REAL :: tem2d7 (m,n)
  REAL :: tem2d8 (m,n)
  REAL :: tem2d9 (m,n)
  REAL :: tem2d10(m,n)
  REAL :: tem2d11(m,n)
  REAL :: tem2d12(m,n)
  REAL :: tem2d13(m,n)
  REAL :: tem2d14(m,n)
  REAL :: tem2d15(m,n)
  REAL :: tem2d16(m,n)
  REAL :: tem2d17(m,n)
  REAL :: tem2d18(m,n)
  REAL :: tem2d19(m,n)

  REAL :: tem3d1(m,n,np)
  REAL :: tem3d2(m,n,np)
  REAL :: tem3d3(m,n,np+1)
  REAL :: tem3d4(m,n,np+1)
  REAL :: tem3d5(m,n,np+1)

  REAL :: tem4d1(m,n,np+1,2)
  REAL :: tem4d2(m,n,np+1,2)
  REAL :: tem4d3(m,n,np+1,2)
  REAL :: tem4d4(m,n,np+1,2)
  REAL :: tem4d5(m,n,np+1,2)

  REAL :: tem5d1(m,n,np+1,2,2)
  REAL :: tem5d2(m,n,np+1,2,2)
  REAL :: tem5d3(m,n,np+1,2,2)
  REAL :: tem5d4(m,n,np+1,2,2)
  REAL :: tem5d5(m,n,np+1,2,2)

!-----temporary variables

  INTEGER :: i,j,k
  REAL :: x
!
!-----------------------------------------------------------------------
!
!  Functions:
!
!-----------------------------------------------------------------------
!
  REAL :: expmn
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO j=1,n
    DO i=1,m

      swh(i,j,1)=0.
      so2(i,j,1)=0.

!-----csm is the effective secant of the solar zenith angle
!  see equation (12) of Lacis and Hansen (1974, JAS)

      csm(i,j)=35./SQRT(1224.*cosz(i,j)*cosz(i,j)+1.)

    END DO
  END DO


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

!-----compute layer thickness and pressure-scaling function.
!  indices for the surface level and surface layer
!  are np+1 and np, respectively.

        dp(i,j,k)=pl(i,j,k+1)-pl(i,j,k)
        scal(i,j,k)=dp(i,j,k)*(.5*(pl(i,j,k)+pl(i,j,k+1))/300.)**.8

!-----compute scaled water vapor amount, unit is g/cm**2

        wh(i,j,k)=1.02*wa(i,j,k)*scal(i,j,k)*                           &
                  (1.-0.00135*(ta(i,j,k)-240.))
        swh(i,j,k+1)=swh(i,j,k)+wh(i,j,k)

!-----compute ozone amount, unit is (cm-atm)stp.

        oh(i,j,k)=1.02*oa(i,j,k)*dp(i,j,k)*466.7

      END DO
    END DO
  END DO


!-----scale cloud optical thickness in each layer from taucld (with
!  cloud amount fcld) to tauclb and tauclf (with cloud amount cc).
!  tauclb is the scaled optical thickness for beam radiation and
!  tauclf is for diffuse radiation.

  CALL cldscale(idim,jdim,                                              &
                m,n,np,cosz,fcld,taucldi,taucldl,ict,icb,               &
                cc,tauclb,tauclf)

!-----initialize fluxes for all-sky (flx), clear-sky (flc), and
!  flux reduction (df)

  DO k=1,np+1
    DO j=1,n
      DO i=1,m
        flx(i,j,k)=0.
        flc(i,j,k)=0.
        df(i,j,k)=0.
      END DO
    END DO
  END DO

!-----compute solar ir fluxes

  CALL solir (idim,jdim,                                                &
              m,n,np,wh,taucldi,taucldl,tauclb,tauclf,                  &
              reffi,reffl,ict,icb,                                      &
              fcld,cc,taual,csm,rsirbm,rsirdf,                          &
              flx,flc,fdirir,fdifir,                                    &
              tem2d1, tem2d2, tem2d3, tem2d4,                           &
              tem2d5, tem2d6, tem2d7, tem2d8, tem2d9,                   &
              tem2d10,tem2d11,tem2d12,tem2d13,tem2d14,                  &
              tem3d1, tem3d2, tem3d3, tem3d4,                           &
              tem4d1, tem4d2, tem4d3, tem4d4, tem4d5,                   &
              tem2d15,tem2d16,tem2d17,tem2d18,tem2d19,                  &
              tem3d5,                                                   &
              tem5d1,tem5d2,tem5d3,tem5d4,tem5d5)

!  DO k=1,np
!    write (6,'(a,i2,a,i2,a,e20.10)')
!    :  'IR:  flx(2,2,',k+1,')-flx(2,2,',k,')=',
!    :  flx(2,2,k+1)-flx(2,2,k)
!  ENDDO

!-----compute and update uv and par fluxes

  CALL soluv (idim,jdim,                                                &
              m,n,np,oh,dp,taucldi,taucldl,tauclb,tauclf,               &
              reffi,reffl,ict,icb,                                      &
              fcld,cc,taual,csm,rsuvbm,rsuvdf,                          &
              flx,flc,fdirpar,fdifpar,                                  &
              tem2d1, tem2d2, tem2d3,                                   &
              tem2d5, tem2d6, tem2d7, tem2d8, tem2d9,                   &
              tem2d10,tem2d11,tem2d12,tem2d13,tem2d14,                  &
              tem3d1, tem3d3, tem3d4,                                   &
              tem4d1, tem4d2, tem4d3, tem4d4, tem4d5,                   &
              tem2d15,tem2d16,tem2d17,tem2d18,tem2d19,                  &
              tem3d5,                                                   &
              tem5d1,tem5d2,tem5d3,tem5d4,tem5d5)

!  write (6,'(a/a3,2a15)') 'Total flux for PAR',
!    :  'k','dflxPAR','flxPAR'
!  DO k=1,np
!    write (6,'(i3,2e15.7)')
!    :  k,flx(2,2,k+1)-flx(2,2,k),flx(2,2,k)
!  ENDDO

!-----compute scaled amount of o2 (so2), unit is (cm-atm)stp.

  DO k=1,np
    DO j=1,n
      DO i=1,m
        so2(i,j,k+1)=so2(i,j,k)+165.22*scal(i,j,k)
      END DO
    END DO
  END DO

!-----compute flux reduction due to oxygen following
!   chou (J. climate, 1990). The fraction 0.0287 is the
!   extraterrestrial solar flux in the o2 bands.

  DO k=2,np+1
    DO j=1,n
      DO i=1,m
        x=so2(i,j,k)*csm(i,j)
        df(i,j,k)=df(i,j,k)+0.0287*(1.-expmn(-0.00027*SQRT(x)))
      END DO
    END DO
  END DO

!  write (6,'(a/a3,3a15)')
!    : 'Flux reduction due to oxygen, i=j=3',
!    : 'k','ddf','df','so2'
!  DO k=2,np+1
!    write (6,'(i3,3e15.7)')
!    :  k,df(3,3,k)-df(3,3,k-1),df(3,3,k),so2(3,3,k)
!  ENDDO

!-----compute scaled amounts for co2 (so2). unit is (cm-atm)stp.

  DO k=1,np
    DO j=1,n
      DO i=1,m
        so2(i,j,k+1)=so2(i,j,k)+co2*789.*scal(i,j,k)
      END DO
    END DO
  END DO

!-----compute and update flux reduction due to co2 following
!  chou (J. Climate, 1990)

  CALL flxco2(m,n,np,so2,swh,csm,df)

!-----adjust for the effect of o2 cnd co2 on clear-sky fluxes.

  DO k=2,np+1
    DO j=1,n
      DO i=1,m
        flc(i,j,k)=flc(i,j,k)-df(i,j,k)
      END DO
    END DO
  END DO

!-----adjust for the all-sky fluxes due to o2 and co2.  It is
!  assumed that o2 and co2 have no effects on solar radiation
!  below clouds. clouds are assumed randomly overlapped.

  DO j=1,n
    DO i=1,m
      sdf(i,j)=0.0
      sclr(i,j)=1.0
    END DO
  END DO

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

!-----sclr is the fraction of clear sky.
!  sdf is the flux reduction below clouds.

        IF(fcld(i,j,k) > 0.01) THEN
          sdf(i,j)=sdf(i,j)+df(i,j,k)*sclr(i,j)*fcld(i,j,k)
          sclr(i,j)=sclr(i,j)*(1.-fcld(i,j,k))
        END IF
        flx(i,j,k+1)=flx(i,j,k+1)-sdf(i,j)-df(i,j,k+1)*sclr(i,j)

      END DO
    END DO
  END DO

!-----adjust for the direct downward ir flux.
  DO j=1,n
    DO i=1,m
      fdirir(i,j)=fdirir(i,j)-sdf(i,j)-df(i,j,np+1)*sclr(i,j)
    END DO
  END DO

!  write (6,'(a3,2a15)') 'k','flx','flc'
!  DO k=1,np
!    write (6,'(i3,2e15.8)') k,flx(1,1,k),flc(1,1,k)
!  ENDDO
!
!  write (6,'(4a15)')
!    : 'fdirir','fdifir','fdirpar','fdifpar'
!
!  write (6,'(4e15.8)')
!    :  fdirir(1,1),fdifir(1,1),fdirpar(1,1),fdifpar(1,1)

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


SUBROUTINE solir(idim,jdim,                                             & 1,5
           m,n,np,wh,taucldi,taucldl,tauclb,tauclf,reffi,reffl,         &
           ict,icb,fcld,cc,taual,csm,rsirbm,rsirdf,                     &
           flx,flc,fdirir,fdifir,                                       &
           fsdir,fsdif, ssaclt,asyclt,                                  &
           rr1t,tt1t,td1t,rs1t,ts1t,                                    &
           rr2t,tt2t,td2t,rs2t,ts2t,                                    &
           ssacl,asycl,fall,fclr,                                       &
           rr,tt,td,rs,ts,                                              &
           tem2d1,tem2d2,tem2d3,tem2d4,tem2d5,                          &
           tem3d,                                                       &
           tem5d1,tem5d2,tem5d3,tem5d4,tem5d5)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate solar flux in the infrared region. The spectrum is
!  divided into three bands. (See original comments.)
!
!-----------------------------------------------------------------------
!
!  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)
!  Modified the original code from 1-D to 3-D
!
!  12/14/1998 (Keith Brewster)
!  Modified values of ssaal and asymal in data statement.
!
!-----------------------------------------------------------------------
!
!  ORIGINAL COMMENTS:
!
!************************************************************************
!  compute solar flux in the infrared region. The spectrum is divided
!   into three bands:
!
!       band   wavenumber(/cm)  wavelength (micron)
!        1       1000-4400         2.27-10.0
!        2       4400-8200         1.22-2.27
!        3       8200-14300        0.70-1.22
!
!----- 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                   n/d        1
!       meridional direction (ndim)
!   number of atmospheric layers (np)                n/d        1
!   layer water vapor content (wh)                 gm/cm^2    m*n*np
!   cloud optical thickness (taucld)                 n/d      m*ndim*np*2
!       index 1 for ice paticles
!       index 2 for liquid particles
!   scaled cloud optical thickness                   n/d      m*n*np
!       for beam radiation (tauclb)
!   scaled cloud optical thickness                   n/d      m*n*np
!       for diffuse radiation  (tauclf)
!   effective cloud-particle size (reff)           micrometer m*ndim*np*2
!       index 1 for ice paticles
!       index 2 for liquid particles
!   level index separating high and                  n/d      m*n
!       middle clouds (ict)
!   level index separating middle and                n/d      m*n
!       low clouds (icb)
!   cloud amount (fcld)                            fraction   m*ndim*np
!   cloud amount of high, middle, and                n/d      m*n*3
!       low clouds (cc)
!   aerosol optical thickness (taual)                n/d      m*ndim*np
!   cosecant of the solar zenith angle (csm)         n/d      m*n
!   near ir surface albedo for beam                fraction   m*ndim
!             radiation (rsirbm)
!   near ir surface albedo for diffuse             fraction   m*ndim
!             radiation (rsirdf)
!
!----- output (updated) parameters:
!
!   all-sky flux (downward-upward) (flx)           fraction   m*ndim*(np+1)
!   clear-sky flux (downward-upward) (flc)         fraction   m*ndim*(np+1)
!   all-sky direct downward ir flux at
!       the surface (fdirir)                    fraction   m*ndim
!   all-sky diffuse downward ir flux at
!       the surface (fdifir)                    fraction   m*ndim
!
!----- note: the following parameters must be specified by users:
!   aerosol single scattering albedo (ssaal)         n/d      nband
!   aerosol asymmetry factor (asyal)                 n/d      nband
!
!*************************************************************************
!
!fpp$ expand (expmn)
!fpp$ expand (deledd)
!fpp$ expand (sagpol)
!dir$ inline always expmn,deledd,sagpol
!*$*  inline routine (expmn,deledd,sagpol)
!
!-----------------------------------------------------------------------
!
!  Variable declarations
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

!-----input parameters

  INTEGER :: idim,jdim
  INTEGER :: m,n,np,ict,icb

  REAL :: taucldi(idim,jdim,np)
  REAL :: taucldl(idim,jdim,np)
  REAL :: reffi(idim,jdim,np)
  REAL :: reffl(idim,jdim,np)
  REAL :: fcld(idim,jdim,np)
  REAL :: rsirbm(idim,jdim)
  REAL :: rsirdf(idim,jdim)
  REAL :: taual(idim,jdim,np)

  REAL :: tauclb(m,n,np)
  REAL :: tauclf(m,n,np)
  REAL :: cc(m,n,3)
  REAL :: wh(m,n,np)
  REAL :: csm(m,n)

!-----output (updated) parameters

  REAL :: flx(idim,jdim,np+1)
  REAL :: flc(idim,jdim,np+1)
  REAL :: fdirir(idim,jdim)
  REAL :: fdifir(idim,jdim)

!-----static parameters

  INTEGER :: nk,nband
  PARAMETER (nk=10,nband=3)
  REAL :: xk(nk),hk(nband,nk),ssaal(nband),asyal(nband)
  REAL :: aia(nband,3),awa(nband,3),aig(nband,3),awg(nband,3)

!-----temporary array

  REAL :: fsdir(m,n)
  REAL :: fsdif(m,n)

  REAL :: ssaclt(m,n)
  REAL :: asyclt(m,n)
  REAL :: rr1t(m,n)
  REAL :: tt1t(m,n)
  REAL :: td1t(m,n)
  REAL :: rs1t(m,n)
  REAL :: ts1t(m,n)
  REAL :: rr2t(m,n)
  REAL :: tt2t(m,n)
  REAL :: td2t(m,n)
  REAL :: rs2t(m,n)
  REAL :: ts2t(m,n)

  REAL :: ssacl(m,n,np)
  REAL :: asycl(m,n,np)
  REAL :: fall(m,n,np+1)
  REAL :: fclr(m,n,np+1)

  REAL :: rr(m,n,np+1,2)
  REAL :: tt(m,n,np+1,2)
  REAL :: td(m,n,np+1,2)
  REAL :: rs(m,n,np+1,2)
  REAL :: ts(m,n,np+1,2)

!-----temporary arrays used in subroutine cldflx

  REAL :: tem2d1(m,n)
  REAL :: tem2d2(m,n)
  REAL :: tem2d3(m,n)
  REAL :: tem2d4(m,n)
  REAL :: tem2d5(m,n)

  REAL :: tem3d (m,n,np+1)

  REAL :: tem5d1(m,n,np+1,2,2)
  REAL :: tem5d2(m,n,np+1,2,2)
  REAL :: tem5d3(m,n,np+1,2,2)
  REAL :: tem5d4(m,n,np+1,2,2)
  REAL :: tem5d5(m,n,np+1,2,2)

!-----temporary variables

  INTEGER :: ib,ik,i,j,k

  REAL :: tauwv,tausto,ssatau,asysto,tauto,ssato,asyto
  REAL :: taux,reff1,reff2,w1,w2,g1,g2

!-----function

  REAL :: expmn

!-----water vapor absorption coefficient for 10 k-intervals.
!  unit: cm^2/gm

  DATA xk/                                                              &
      0.0010, 0.0133, 0.0422, 0.1334, 0.4217,                           &
      1.334,  5.623,  31.62,  177.8,  1000.0/

!-----water vapor k-distribution function,
!  the sum of hk is 0.52926. unit: fraction

  DATA hk/                                                              &
      .01074,.08236,.20673,  .00360,.01157,.03497,                      &
      .00411,.01133,.03011,  .00421,.01143,.02260,                      &
      .00389,.01240,.01336,  .00326,.01258,.00696,                      &
      .00499,.01381,.00441,  .00465,.00650,.00115,                      &
      .00245,.00244,.00026,  .00145,.00094,.00000/

!-----aerosol single-scattering albedo and asymmetry factor

!  data ssaal,asyal/nband*0.999,nband*0.850/
  DATA ssaal/0.75,0.55,0.90/
  DATA asyal/0.40,0.50,0.60/

!  data ssaal,asyal/nband*0.400,nband*0.500/

!-----coefficients for computing the single scattering albedo of
!  ice clouds from ssa=1-(aia(*,1)+aia(*,2)*reff+aia(*,3)*reff**2)

  DATA aia/                                                             &
      .08938331, .00215346,-.00000260,                                  &
      .00299387, .00073709, .00000746,                                  &
      -.00001038,-.00000134, .00000000/

!-----coefficients for computing the single scattering albedo of
!  liquid clouds from ssa=1-(awa(*,1)+awa(*,2)*reff+awa(*,3)*reff**2)

  DATA awa/                                                             &
      .01209318,-.00019934, .00000007,                                  &
      .01784739, .00088757, .00000845,                                  &
      -.00036910,-.00000650,-.00000004/

!-----coefficients for computing the asymmetry factor of ice clouds
!  from asycl=aig(*,1)+aig(*,2)*reff+aig(*,3)*reff**2

  DATA aig/                                                             &
      .84090400, .76098937, .74935228,                                  &
      .00126222, .00141864, .00119715,                                  &
      -.00000385,-.00000396,-.00000367/

!-----coefficients for computing the asymmetry factor of liquid clouds
!  from asycl=awg(*,1)+awg(*,2)*reff+awg(*,3)*reff**2

  DATA awg/                                                             &
      .83530748, .74513197, .79375035,                                  &
      .00257181, .01370071, .00832441,                                  &
      .00005519,-.00038203,-.00023263/
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

!-----initialize surface fluxes, reflectances, and transmittances

  DO j= 1, n
    DO i= 1, m
      fdirir(i,j)=0.0
      fdifir(i,j)=0.0
      rr(i,j,np+1,1)=rsirbm(i,j)
      rr(i,j,np+1,2)=rsirdf(i,j)
      rs(i,j,np+1,1)=rsirbm(i,j)
      rs(i,j,np+1,2)=rsirdf(i,j)
      td(i,j,np+1,1)=0.0
      td(i,j,np+1,2)=0.0
      tt(i,j,np+1,1)=0.0
      tt(i,j,np+1,2)=0.0
      ts(i,j,np+1,1)=0.0
      ts(i,j,np+1,2)=0.0
    END DO
  END DO

!-----integration over spectral bands

  DO ib=1,nband

!-----compute cloud single scattering albedo and asymmetry factor
!  for a mixture of ice and liquid particles.

    DO k= 1, np

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

          ssaclt(i,j)=1.0
          asyclt(i,j)=1.0

          taux=taucldi(i,j,k)+taucldl(i,j,k)
          IF (taux > 0.05 .AND. fcld(i,j,k) > 0.01) THEN

            reff1=MIN(reffi(i,j,k),130.)
            reff2=MIN(reffl(i,j,k),20.0)

            w1=(1.-(aia(ib,1)+(aia(ib,2)+                               &
                aia(ib,3)*reff1)*reff1))*taucldi(i,j,k)
            w2=(1.-(awa(ib,1)+(awa(ib,2)+                               &
                awa(ib,3)*reff2)*reff2))*taucldl(i,j,k)
            ssaclt(i,j)=(w1+w2)/taux

            g1=(aig(ib,1)+(aig(ib,2)+aig(ib,3)*reff1)*reff1)*w1
            g2=(awg(ib,1)+(awg(ib,2)+awg(ib,3)*reff2)*reff2)*w2
            asyclt(i,j)=(g1+g2)/(w1+w2)

          END IF

        END DO
      END DO

      DO j=1,n
        DO i=1,m
          ssacl(i,j,k)=ssaclt(i,j)
        END DO
      END DO
      DO j=1,n
        DO i=1,m
          asycl(i,j,k)=asyclt(i,j)
        END DO
      END DO

    END DO

!-----integration over the k-distribution function

    DO ik=1,nk

      DO k= 1, np

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

            tauwv=xk(ik)*wh(i,j,k)

!-----compute total optical thickness, single scattering albedo,
!  and asymmetry factor.

            tausto=tauwv+taual(i,j,k)+1.0E-8
            ssatau=ssaal(ib)*taual(i,j,k)
            asysto=asyal(ib)*ssaal(ib)*taual(i,j,k)

!-----compute reflectance and transmittance for cloudless layers

            tauto=tausto
            ssato=ssatau/tauto+1.0E-8

            IF (ssato > 0.001) THEN

              ssato=MIN(ssato,0.999999)
              asyto=asysto/(ssato*tauto)

              CALL deledd (tauto,ssato,asyto,csm(i,j),                  &
                           rr1t(i,j),tt1t(i,j),td1t(i,j))

              CALL sagpol (tauto,ssato,asyto,rs1t(i,j),ts1t(i,j))

            ELSE

              td1t(i,j)=expmn(-tauto*csm(i,j))
              ts1t(i,j)=expmn(-1.66*tauto)
              tt1t(i,j)=0.0
              rr1t(i,j)=0.0
              rs1t(i,j)=0.0

            END IF

!-----compute reflectance and transmittance for cloud layers

            IF (tauclb(i,j,k) < 0.01) THEN

              rr2t(i,j)=rr1t(i,j)
              tt2t(i,j)=tt1t(i,j)
              td2t(i,j)=td1t(i,j)
              rs2t(i,j)=rs1t(i,j)
              ts2t(i,j)=ts1t(i,j)

            ELSE

              tauto=tausto+tauclb(i,j,k)
              ssato=(ssatau+ssacl(i,j,k)*tauclb(i,j,k))/tauto+1.0E-8
              ssato=MIN(ssato,0.999999)
              asyto=(asysto+asycl(i,j,k)*ssacl(i,j,k)*tauclb(i,j,k))/   &
                    (ssato*tauto)

              CALL deledd (tauto,ssato,asyto,csm(i,j),                  &
                           rr2t(i,j),tt2t(i,j),td2t(i,j))

              tauto=tausto+tauclf(i,j,k)
              ssato=(ssatau+ssacl(i,j,k)*tauclf(i,j,k))/tauto+1.0E-8
              ssato=MIN(ssato,0.999999)
              asyto=(asysto+asycl(i,j,k)*ssacl(i,j,k)*tauclf(i,j,k))/   &
                    (ssato*tauto)

              CALL sagpol (tauto,ssato,asyto,rs2t(i,j),ts2t(i,j))

            END IF

          END DO
        END DO

        DO j=1,n
          DO i=1,m
            rr(i,j,k,1)=rr1t(i,j)
          END DO
        END DO
        DO j=1,n
          DO i=1,m
            tt(i,j,k,1)=tt1t(i,j)
          END DO
        END DO
        DO j=1,n
          DO i=1,m
            td(i,j,k,1)=td1t(i,j)
          END DO
        END DO
        DO j=1,n
          DO i=1,m
            rs(i,j,k,1)=rs1t(i,j)
          END DO
        END DO
        DO j=1,n
          DO i=1,m
            ts(i,j,k,1)=ts1t(i,j)
          END DO
        END DO

        DO j=1,n
          DO i=1,m
            rr(i,j,k,2)=rr2t(i,j)
          END DO
        END DO
        DO j=1,n
          DO i=1,m
            tt(i,j,k,2)=tt2t(i,j)
          END DO
        END DO
        DO j=1,n
          DO i=1,m
            td(i,j,k,2)=td2t(i,j)
          END DO
        END DO
        DO j=1,n
          DO i=1,m
            rs(i,j,k,2)=rs2t(i,j)
          END DO
        END DO
        DO j=1,n
          DO i=1,m
            ts(i,j,k,2)=ts2t(i,j)
          END DO
        END DO

      END DO

!-----flux calculations

      CALL cldflx (m,n,np,ict,icb,cc,rr,tt,td,rs,ts,                    &
                   fclr,fall,fsdir,fsdif,                               &
                   tem2d1,tem2d2,tem2d3,tem2d4,tem2d5,                  &
                   tem3d, tem5d1,tem5d2,tem5d3,tem5d4,tem5d5)

      DO k= 1, np+1
        DO j= 1, n
          DO i= 1, m
            flx(i,j,k) = flx(i,j,k)+fall(i,j,k)*hk(ib,ik)
          END DO
        END DO
        DO j= 1, n
          DO i= 1, m
            flc(i,j,k) = flc(i,j,k)+fclr(i,j,k)*hk(ib,ik)
          END DO
        END DO
      END DO

!  write (6,'(a/2a3,2a15)') 'Flux for IR',
!    :  'ib','k','dflxIR','flxIR'
!  DO k=1,np
!    write (6,'(2i3,2e15.7)')
!    :  ib,k,fall(2,2,k+1)-fall(2,2,k),fall(2,2,k)
!  ENDDO

      DO j= 1, n
        DO i= 1, m
          fdirir(i,j) = fdirir(i,j)+fsdir(i,j)*hk(ib,ik)
          fdifir(i,j) = fdifir(i,j)+fsdif(i,j)*hk(ib,ik)
        END DO
      END DO

    END DO
  END DO

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


SUBROUTINE soluv(idim,jdim,                                             & 1,5
           m,n,np,oh,dp,taucldi,taucldl,tauclb,tauclf,reffi,reffl,      &
           ict,icb,fcld,cc,taual,csm,rsuvbm,rsuvdf,                     &
           flx,flc,fdirpar,fdifpar,                                     &
           fsdir,fsdif,asyclt,                                          &
           rr1t,tt1t,td1t,rs1t,ts1t,                                    &
           rr2t,tt2t,td2t,rs2t,ts2t,                                    &
           asycl,fall,fclr,                                             &
           td,rr,tt,rs,ts,                                              &
           tem2d1,tem2d2,tem2d3,tem2d4,tem2d5,                          &
           tem3d,                                                       &
           tem5d1,tem5d2,tem5d3,tem5d4,tem5d5)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate solar fluxes in the uv+visible region. the spectrum is
!  grouped into 8 bands. (See original comments.)
!
!-----------------------------------------------------------------------
!
!  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)
!  Adapted the original code and formatted it in accordance with the
!  ARPS coding standard.
!
!  12/07/1998 (Keith Brewster)
!  Modified values of ssaal and asymal in data statement.
!
!-----------------------------------------------------------------------
!
!  ORIGINAL COMMENTS:
!
!************************************************************************
!  compute solar fluxes in the uv+visible region. the spectrum is
!  grouped into 8 bands:
!
!           Band     Micrometer
!
!    UV-C    1.     .175 - .225
!            2.     .225 - .245
!                   .260 - .280
!            3.     .245 - .260
!
!    UV-B    4.     .280 - .295
!            5.     .295 - .310
!            6.     .310 - .320
!
!    UV-A    7.     .320 - .400
!
!    PAR     8.     .400 - .700
!
!----- 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                   n/d        1
!        meridional direction (ndim)
!   number of atmospheric layers (np)                n/d        1
!   layer ozone content (oh)                      (cm-atm)stp m*n*np
!   layer pressure thickness (dp)                    mb       m*n*np
!   cloud optical thickness (taucld)                 n/d      m*ndim*np*2
!       index 1 for ice paticles
!       index 2 for liquid particles
!   scaled cloud optical thickness                   n/d      m*n*np
!       for beam radiation (tauclb)
!   scaled cloud optical thickness                   n/d      m*n*np
!       for diffuse radiation  (tauclf)
!   effective cloud-particle size (reff)           micrometer m*ndim*np*2
!       index 1 for ice paticles
!       index 2 for liquid particles
!   level indiex separating high and                 n/d      m*n
!       middle clouds (ict)
!   level indiex separating middle and               n/d      m*n
!       low clouds (icb)
!   cloud amount (fcld)                            fraction   m*ndim*np
!   cloud amounts of high, middle, and               n/d      m*n*3
!       low clouds (cc)
!   aerosol optical thickness (taual)                n/d      m*ndim*np
!   cosecant of the solar zenith angle (csm)         n/d      m*n
!   uv+par surface albedo for beam                 fraction   m*ndim
!        radiation (rsuvbm)
!   uv+par surface albedo for diffuse              fraction   m*ndim
!        radiation (rsuvdf)
!
!----- output (updated) parameters:
!
!   all-sky net downward flux (flx)                fraction   m*ndim*(np+1)
!   clear-sky net downward flux (flc)              fraction   m*ndim*(np+1)
!   all-sky direct downward par flux at
!       the surface (fdirpar)                   fraction   m*ndim
!   all-sky diffuse downward par flux at
!       the surface (fdifpar)                   fraction   m*ndim
!
!----- note: the following parameters must be specified by users:
!
!   aerosol single scattering albedo (ssaal)         n/d        1
!   aerosol asymmetry factor (asyal)                 n/d        1
!
!
!***********************************************************************
!
!fpp$ expand (deledd)
!fpp$ expand (sagpol)
!dir$ inline always deledd,sagpol
!*$*  inline routine (deledd,sagpol)
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

!-----input parameters

  INTEGER :: idim,jdim
  INTEGER :: m,n,np,ict,icb

  REAL :: taucldi(idim,jdim,np)
  REAL :: taucldl(idim,jdim,np)
  REAL :: reffi(idim,jdim,np)
  REAL :: reffl(idim,jdim,np)
  REAL :: fcld(idim,jdim,np)
  REAL :: taual(idim,jdim,np)
  REAL :: rsuvbm(idim,jdim)
  REAL :: rsuvdf(idim,jdim)

  REAL :: tauclb(m,n,np)
  REAL :: tauclf(m,n,np)
  REAL :: cc(m,n,3)
  REAL :: oh(m,n,np)
  REAL :: dp(m,n,np)
  REAL :: csm(m,n)

!-----output (updated) parameter

  REAL :: flx(idim,jdim,np+1)
  REAL :: flc(idim,jdim,np+1)
  REAL :: fdirpar(idim,jdim)
  REAL :: fdifpar(idim,jdim)

!-----static parameters

  INTEGER :: nband
  PARAMETER (nband=8)
  REAL :: hk(nband)
  REAL :: xk(nband)
  REAL :: ry(nband)
  REAL :: asyal(nband)
  REAL :: ssaal(nband)
  REAL :: aig(3)
  REAL :: awg(3)

!-----temporary array

  REAL :: fsdir(m,n)
  REAL :: fsdif(m,n)
  REAL :: asyclt(m,n)
  REAL :: rr1t(m,n)
  REAL :: tt1t(m,n)
  REAL :: td1t(m,n)
  REAL :: rs1t(m,n)
  REAL :: ts1t(m,n)
  REAL :: rr2t(m,n)
  REAL :: tt2t(m,n)
  REAL :: td2t(m,n)
  REAL :: rs2t(m,n)
  REAL :: ts2t(m,n)

  REAL :: asycl(m,n,np)
  REAL :: fall(m,n,np+1)
  REAL :: fclr(m,n,np+1)

  REAL :: td(m,n,np+1,2)
  REAL :: rr(m,n,np+1,2)
  REAL :: tt(m,n,np+1,2)
  REAL :: rs(m,n,np+1,2)
  REAL :: ts(m,n,np+1,2)

!-----temporary arrays used in subroutine cldflx

  REAL :: tem2d1(m,n)
  REAL :: tem2d2(m,n)
  REAL :: tem2d3(m,n)
  REAL :: tem2d4(m,n)
  REAL :: tem2d5(m,n)

  REAL :: tem3d (m,n,np+1)

  REAL :: tem5d1(m,n,np+1,2,2)
  REAL :: tem5d2(m,n,np+1,2,2)
  REAL :: tem5d3(m,n,np+1,2,2)
  REAL :: tem5d4(m,n,np+1,2,2)
  REAL :: tem5d5(m,n,np+1,2,2)

!-----temporary variables

  INTEGER :: i,j,k,ib
  REAL :: taurs,tauoz,tausto,ssatau,asysto,tauto,ssato,asyto
  REAL :: taux,reff1,reff2,g1,g2

!-----hk is the fractional extra-terrestrial solar flux.
!  the sum of hk is 0.47074.

  DATA hk/.00057, .00367, .00083, .00417,                               &
          .00600, .00556, .05913, .39081/

!-----xk is the ozone absorption coefficient. unit: /(cm-atm)stp

  DATA xk /30.47, 187.2,  301.9,   42.83,                               &
           7.09,  1.25,   0.0345,  0.0539/

!-----ry is the extinction coefficient for Rayleigh scattering.
!  unit: /mb.

  DATA ry /.00604, .00170, .00222, .00132,                              &
           .00107, .00091, .00055, .00012/

!-----aerosol single-scattering albedo and asymmetry factor

!  data ssaal,asyal/nband*0.999,nband*0.850/
  DATA ssaal,asyal/nband*0.960,nband*0.700/

!-----coefficients for computing the asymmetry factor of ice clouds
!  from asycl=aig(*,1)+aig(*,2)*reff+aig(*,3)*reff**2

  DATA aig/.74625000,.00105410,-.00000264/

!-----coefficients for computing the asymmetry factor of liquid
!  clouds from asycl=awg(*,1)+awg(*,2)*reff+awg(*,3)*reff**2

  DATA awg/.82562000,.00529000,-.00014866/

!-----initialize surface reflectances and transmittances

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO j= 1, n
    DO i= 1, m
      rr(i,j,np+1,1)=rsuvbm(i,j)
      rr(i,j,np+1,2)=rsuvdf(i,j)
      rs(i,j,np+1,1)=rsuvbm(i,j)
      rs(i,j,np+1,2)=rsuvdf(i,j)
      td(i,j,np+1,1)=0.0
      td(i,j,np+1,2)=0.0
      tt(i,j,np+1,1)=0.0
      tt(i,j,np+1,2)=0.0
      ts(i,j,np+1,1)=0.0
      ts(i,j,np+1,2)=0.0
    END DO
  END DO

!-----compute cloud asymmetry factor for a mixture of
!  liquid and ice particles.  unit of reff is micrometers.

  DO k= 1, np

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

        asyclt(i,j)=1.0

        taux=taucldi(i,j,k)+taucldl(i,j,k)
        IF (taux > 0.05 .AND. fcld(i,j,k) > 0.01) THEN

          reff1=MIN(reffi(i,j,k),130.)
          reff2=MIN(reffl(i,j,k),20.0)

          g1=(aig(1)+(aig(2)+aig(3)*reff1)*reff1)*taucldi(i,j,k)
          g2=(awg(1)+(awg(2)+awg(3)*reff2)*reff2)*taucldl(i,j,k)
          asyclt(i,j)=(g1+g2)/taux

        END IF

      END DO
    END DO

    DO j=1,n
      DO i=1,m
        asycl(i,j,k)=asyclt(i,j)
      END DO
    END DO

  END DO

!-----integration over spectral bands

  DO ib=1,nband

    DO k= 1, np

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

!-----compute ozone and rayleigh optical thicknesses

          taurs=ry(ib)*dp(i,j,k)
          tauoz=xk(ib)*oh(i,j,k)

!-----compute total optical thickness, single scattering albedo,
!  and asymmetry factor

          tausto=taurs+tauoz+taual(i,j,k)+1.0E-8
          ssatau=ssaal(ib)*taual(i,j,k)+taurs
          asysto=asyal(ib)*ssaal(ib)*taual(i,j,k)

!-----compute reflectance and transmittance for cloudless layers

          tauto=tausto
          ssato=ssatau/tauto+1.0E-8
          ssato=MIN(ssato,0.999999)
          asyto=asysto/(ssato*tauto)

          CALL deledd (tauto,ssato,asyto,csm(i,j),                      &
                       rr1t(i,j),tt1t(i,j),td1t(i,j))

          CALL sagpol (tauto,ssato,asyto,rs1t(i,j),ts1t(i,j))

!-----compute reflectance and transmittance for cloud layers

          IF (tauclb(i,j,k) < 0.01) THEN

            rr2t(i,j)=rr1t(i,j)
            tt2t(i,j)=tt1t(i,j)
            td2t(i,j)=td1t(i,j)
            rs2t(i,j)=rs1t(i,j)
            ts2t(i,j)=ts1t(i,j)

          ELSE

            tauto=tausto+tauclb(i,j,k)
            ssato=(ssatau+tauclb(i,j,k))/tauto+1.0E-8
            ssato=MIN(ssato,0.999999)
            asyto=(asysto+asycl(i,j,k)*tauclb(i,j,k))/(ssato*tauto)

            CALL deledd (tauto,ssato,asyto,csm(i,j),                    &
                         rr2t(i,j),tt2t(i,j),td2t(i,j))

            tauto=tausto+tauclf(i,j,k)
            ssato=(ssatau+tauclf(i,j,k))/tauto+1.0E-8
            ssato=MIN(ssato,0.999999)
            asyto=(asysto+asycl(i,j,k)*tauclf(i,j,k))/(ssato*tauto)

            CALL sagpol (tauto,ssato,asyto,rs2t(i,j),ts2t(i,j))

          END IF

        END DO
      END DO

      DO j=1,n
        DO i=1,m
          rr(i,j,k,1)=rr1t(i,j)
        END DO
      END DO
      DO j=1,n
        DO i=1,m
          tt(i,j,k,1)=tt1t(i,j)
        END DO
      END DO
      DO j=1,n
        DO i=1,m
          td(i,j,k,1)=td1t(i,j)
        END DO
      END DO
      DO j=1,n
        DO i=1,m
          rs(i,j,k,1)=rs1t(i,j)
        END DO
      END DO
      DO j=1,n
        DO i=1,m
          ts(i,j,k,1)=ts1t(i,j)
        END DO
      END DO

      DO j=1,n
        DO i=1,m
          rr(i,j,k,2)=rr2t(i,j)
        END DO
      END DO
      DO j=1,n
        DO i=1,m
          tt(i,j,k,2)=tt2t(i,j)
        END DO
      END DO
      DO j=1,n
        DO i=1,m
          td(i,j,k,2)=td2t(i,j)
        END DO
      END DO
      DO j=1,n
        DO i=1,m
          rs(i,j,k,2)=rs2t(i,j)
        END DO
      END DO
      DO j=1,n
        DO i=1,m
          ts(i,j,k,2)=ts2t(i,j)
        END DO
      END DO

    END DO

!-----flux calculations

    CALL cldflx (m,n,np,ict,icb,cc,rr,tt,td,rs,ts,                      &
                 fclr,fall,fsdir,fsdif,                                 &
                 tem2d1,tem2d2,tem2d3,tem2d4,tem2d5,                    &
                 tem3d, tem5d1,tem5d2,tem5d3,tem5d4,tem5d5)

    DO k= 1, np+1
      DO j= 1, n
        DO i= 1, m
          flx(i,j,k)=flx(i,j,k)+fall(i,j,k)*hk(ib)
        END DO
      END DO
      DO j= 1, n
        DO i= 1, m
          flc(i,j,k)=flc(i,j,k)+fclr(i,j,k)*hk(ib)
        END DO
      END DO
    END DO

!  write (6,'(a/2a3,2a15)') 'Flux for PAR',
!    :  'ib','k','dflxpa','flxpa'
!  DO k=1,np
!    write (6,'(2i3,2e15.7)')
!    :  ib,k,fall(2,2,k+1)-fall(2,2,k),fall(2,2,k)
!  ENDDO

    IF(ib == 8) THEN
      DO j=1,n
        DO i=1,m
          fdirpar(i,j)=fsdir(i,j)*hk(ib)
          fdifpar(i,j)=fsdif(i,j)*hk(ib)
        END DO
      END DO
    END IF

  END DO

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


SUBROUTINE cldscale(idim,jdim,                                          & 1
           m,n,np,cosz,fcld,taucldi,taucldl,ict,icb,                    &
           cc,tauclb,tauclf)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate the covers of high, middle, and low clouds and scales
!  the cloud optical thickness.
!
!-----------------------------------------------------------------------
!
!  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:
!
!********************************************************************
!
!   This subroutine computes the covers of high, middle, and
!    low clouds and scales the cloud optical thickness.
!
!   To simplify calculations in a cloudy atmosphere, clouds are
!    grouped into high, middle and low clouds separated by the levels
!    ict and icb (level 1 is the top of the atmosphere).
!
!   Within each of the three groups, clouds are assumed maximally
!    overlapped, and the cloud cover (cc) of a group is the maximum
!    cloud cover of all the layers in the group.  The optical thickness
!    (taucld) of a given layer is then scaled to new values (tauclb and
!    tauclf) so that the layer reflectance corresponding to the cloud
!    cover cc is the same as the original reflectance with optical
!    thickness taucld and cloud cover fcld.
!
!---input parameters
!
!    number of grid intervals in zonal direction (m)
!    number of grid intervals in meridional direction (n)
!    maximum number of grid intervals in meridional direction (ndim)
!    number of atmospheric layers (np)
!    cosine of the solar zenith angle (cosz)
!    fractional cloud cover (fcld)
!    cloud optical thickness (taucld)
!    index separating high and middle clouds (ict)
!    index separating middle and low clouds (icb)
!
!---output parameters
!
!    fractional cover of high, middle, and low clouds (cc)
!    scaled cloud optical thickness for beam radiation (tauclb)
!    scaled cloud optical thickness for diffuse radiation (tauclf)
!
!********************************************************************
!
!-----------------------------------------------------------------------
!
!  Variable declarations
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

!-----input parameters

  INTEGER :: idim,jdim
  INTEGER :: m,n,np,ict,icb

  REAL :: cosz(idim,jdim)
  REAL :: fcld(idim,jdim,np)
  REAL :: taucldi(idim,jdim,np)
  REAL :: taucldl(idim,jdim,np)

!-----output parameters

  REAL :: cc(m,n,3)
  REAL :: tauclb(m,n,np)
  REAL :: tauclf(m,n,np)

!-----temporary variables

  INTEGER :: i,j,k,im,it,ia,kk
  REAL :: fm,ft,fa,xai,taux

!-----pre-computed table

  INTEGER :: nm,nt,na
  PARAMETER (nm=11,nt=9,na=11)
  REAL :: dm,dt,da,t1,caib(nm,nt,na),caif(nt,na)
  PARAMETER (dm=0.1,dt=0.30103,da=0.1,t1=-0.9031)

!-----include the pre-computed table for cai

!  include "cai.dat"
!  save caib,caif

  COMMON /radtab004/ caib,caif
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

!-----clouds within each of the high, middle, and low clouds are
!  assumed maximally overlapped, and the cloud cover (cc)
!  for a group is the maximum cloud cover of all the layers
!  in the group

  DO j=1,n
    DO i=1,m
      cc(i,j,1)=0.0
      cc(i,j,2)=0.0
      cc(i,j,3)=0.0
    END DO
  END DO

  DO k=1,ict-1
    DO j=1,n
      DO i=1,m
        cc(i,j,1)=MAX(cc(i,j,1),fcld(i,j,k))
      END DO
    END DO
  END DO

  DO k=ict,icb-1
    DO j=1,n
      DO i=1,m
        cc(i,j,2)=MAX(cc(i,j,2),fcld(i,j,k))
      END DO
    END DO
  END DO

  DO k=icb,np
    DO j=1,n
      DO i=1,m
        cc(i,j,3)=MAX(cc(i,j,3),fcld(i,j,k))
      END DO
    END DO
  END DO

!-----scale the cloud optical thickness.
!  taucldi(i,j,k) is the optical thickness for ice particles, and
!  taucldl(i,j,k) is the optical thickness for liquid particles.

  DO k=1,np

    IF(k < ict) THEN
      kk=1
    ELSE IF(k >= ict .AND. k < icb) THEN
      kk=2
    ELSE
      kk=3
    END IF

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

        tauclb(i,j,k) = 0.0
        tauclf(i,j,k) = 0.0

        taux=taucldi(i,j,k)+taucldl(i,j,k)
        IF (taux > 0.05 .AND. fcld(i,j,k) > 0.01) THEN

!-----normalize cloud cover

          fa=fcld(i,j,k)/cc(i,j,kk)

!-----table look-up

          taux=MIN(taux,32.)

          fm=cosz(i,j)/dm
          ft=(LOG10(taux)-t1)/dt
          fa=fa/da

          im=INT(fm+1.5)
          it=INT(ft+1.5)
          ia=INT(fa+1.5)

          im=MAX(im,2)
          it=MAX(it,2)
          ia=MAX(ia,2)

          im=MIN(im,nm-1)
          it=MIN(it,nt-1)
          ia=MIN(ia,na-1)

          fm=fm-FLOAT(im-1)
          ft=ft-FLOAT(it-1)
          fa=fa-FLOAT(ia-1)

!-----scale cloud optical thickness for beam radiation.
!  the scaling factor, xai, is a function of the solar zenith
!  angle, optical thickness, and cloud cover.

          xai=    (-caib(im-1,it,ia)*(1.-fm)+                           &
              caib(im+1,it,ia)*(1.+fm))*fm*.5+caib(im,it,ia)*(1.-fm*fm)

          xai=xai+(-caib(im,it-1,ia)*(1.-ft)+                           &
              caib(im,it+1,ia)*(1.+ft))*ft*.5+caib(im,it,ia)*(1.-ft*ft)

          xai=xai+(-caib(im,it,ia-1)*(1.-fa)+                           &
              caib(im,it,ia+1)*(1.+fa))*fa*.5+caib(im,it,ia)*(1.-fa*fa)

          xai= xai-2.*caib(im,it,ia)
          xai=MAX(xai,0.0)

          tauclb(i,j,k) = taux*xai

!-----scale cloud optical thickness for diffuse radiation.
!  the scaling factor, xai, is a function of the cloud optical
!  thickness and cover but not the solar zenith angle.

          xai=    (-caif(it-1,ia)*(1.-ft)+                              &
              caif(it+1,ia)*(1.+ft))*ft*.5+caif(it,ia)*(1.-ft*ft)

          xai=xai+(-caif(it,ia-1)*(1.-fa)+                              &
              caif(it,ia+1)*(1.+fa))*fa*.5+caif(it,ia)*(1.-fa*fa)

          xai= xai-caif(it,ia)
          xai=MAX(xai,0.0)

          tauclf(i,j,k) = taux*xai

        END IF

      END DO
    END DO
  END DO

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


SUBROUTINE flxco2(m,n,np,swc,swh,csm,df) 1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate the reduction of clear-sky downward solar flux
!  due to co2 absorption.
!
!-----------------------------------------------------------------------
!
!  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 the reduction of clear-sky downward solar flux
!  due to co2 absorption.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

!-----input parameters

  INTEGER :: m,n,np

  REAL :: csm(m,n)
  REAL :: swc(m,n,np+1)
  REAL :: swh(m,n,np+1)
  REAL :: cah(22,19)

!-----output (undated) parameter

  REAL :: df(m,n,np+1)

!-----temporary variables

  INTEGER :: i,j,k,ic,iw
  REAL :: xx,CLOG,wlog,dc,dw,x1,x2,y2

!********************************************************************
!-----include co2 look-up table
!
!  include "cah.dat"
!  save cah

  COMMON /radtab005/ cah
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

!********************************************************************
!-----table look-up for the reduction of clear-sky solar
!  radiation due to co2. The fraction 0.0343 is the
!  extraterrestrial solar flux in the co2 bands.

!  write (6,'(a/3a3,7a10)') 'Flux reduction due to co2, i=j=3',
!    :'k','ic','iw','x1','y2','x1-y2','ddf','df','swc','swh'

  DO k= 2, np+1
    DO j= 1, n
      DO i= 1, m
        xx=1./.3
        CLOG=LOG10(swc(i,j,k)*csm(i,j))
        wlog=LOG10(swh(i,j,k)*csm(i,j))
        ic=INT( (CLOG+3.15)*xx+1.)
        iw=INT( (wlog+4.15)*xx+1.)
        IF(ic < 2)ic=2
        IF(iw < 2)iw=2
        IF(ic > 22)ic=22
        IF(iw > 19)iw=19
        dc=CLOG-FLOAT(ic-2)*.3+3.
        dw=wlog-FLOAT(iw-2)*.3+4.
        x1=cah(1,iw-1)+(cah(1,iw)-cah(1,iw-1))*xx*dw
        x2=cah(ic-1,iw-1)+(cah(ic-1,iw)-cah(ic-1,iw-1))*xx*dw
        y2=x2+(cah(ic,iw-1)-cah(ic-1,iw-1))*xx*dc
        IF (x1 < y2) x1=y2
        df(i,j,k)=df(i,j,k)+0.0343*(x1-y2)
      END DO
    END DO

!    write (6,'(3i3,7f10.5)')
!    :  k,ic,iw,x1,y2,x1-y2,
!    :  df(3,3,k)-df(3,3,k-1),df (3,3,k),swc(3,3,k),swh(3,3,k)

  END DO

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


SUBROUTINE cldflx(m,n,np,ict,icb,cc,rr,tt,td,rs,ts,                     & 2
           fclr,fall,fsdir,fsdif,                                       &
           ch,cm,ct,fdndir,fdndif,                                      &
           flxdn,                                                       &
           rra,tta,tda,rsa,rxa)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate upward and downward fluxes using a two-stream adding
!  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 upward and downward fluxes using a two-stream adding method
!  following equations (3)-(5) of Chou (1992, JAS).
!
!  clouds are grouped into high, middle, and low clouds which are
!  assumed randomly overlapped. It involves eight sets of calculations.
!  In each set of calculations, each atmospheric layer is homogeneous,
!  either with clouds or without clouds.

!  input parameters:
!  m:   number of soundings in zonal direction
!  n:   number of soundings in meridional direction
!  np:  number of atmospheric layers
!  ict: the level separating high and middle clouds
!  icb: the level separating middle and low clouds
!  cc:  effective cloud covers for high, middle and low clouds
!  tt:  diffuse transmission of a layer illuminated by beam radiation
!  td:  direct beam tranmssion
!  ts:  transmission of a layer illuminated by diffuse radiation
!  rr:  reflection of a layer illuminated by beam radiation
!  rs:  reflection of a layer illuminated by diffuse radiation
!
!  output parameters:
!  fclr:  clear-sky flux (downward minus upward)
!  fall:  all-sky flux (downward minus upward)
!  fdndir: surface direct downward flux
!  fdndif: surface diffuse downward flux
!*********************************************************************c
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

!-----input parameters

  INTEGER :: m,n,np,ict,icb

  REAL :: rr(m,n,np+1,2)
  REAL :: tt(m,n,np+1,2)
  REAL :: td(m,n,np+1,2)
  REAL :: rs(m,n,np+1,2)
  REAL :: ts(m,n,np+1,2)
  REAL :: cc(m,n,3)

!-----output parameters

  REAL :: fclr(m,n,np+1)
  REAL :: fall(m,n,np+1)
  REAL :: fsdir(m,n)
  REAL :: fsdif(m,n)

!-----temporary array

  REAL :: ch(m,n)
  REAL :: cm(m,n)
  REAL :: ct(m,n)
  REAL :: fdndir(m,n)
  REAL :: fdndif(m,n)

  REAL :: flxdn(m,n,np+1)

  REAL :: rra(m,n,np+1,2,2)
  REAL :: tta(m,n,np+1,2,2)
  REAL :: tda(m,n,np+1,2,2)
  REAL :: rsa(m,n,np+1,2,2)
  REAL :: rxa(m,n,np+1,2,2)

!-----temporary variables

  INTEGER :: i,j,k,ih,im,is
  REAL :: fupdif
  REAL :: denm,xx
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

!-----initialize all-sky flux (fall) and surface downward fluxes

  DO k=1,np+1
    DO j=1,n
      DO i=1,m
        fall(i,j,k)=0.0
      END DO
    END DO
  END DO

  DO j=1,n
    DO i=1,m
      fsdir(i,j)=0.0
      fsdif(i,j)=0.0
    END DO
  END DO

!-----compute transmittances and reflectances for a composite of
!  layers. layers are added one at a time, going down from the top.
!  tda is the composite transmittance illuminated by beam radiation
!  tta is the composite diffuse transmittance illuminated by
!      beam radiation
!  rsa is the composite reflectance illuminated from below
!      by diffuse radiation
!  tta and rsa are computed from eqs. (4b) and (3b) of Chou

!-----for high clouds. indices 1 and 2 denote clear and cloudy
!  situations, respectively.

  DO ih=1,2

    DO j= 1, n
      DO i= 1, m
        tda(i,j,1,ih,1)=td(i,j,1,ih)
        tta(i,j,1,ih,1)=tt(i,j,1,ih)
        rsa(i,j,1,ih,1)=rs(i,j,1,ih)
        tda(i,j,1,ih,2)=td(i,j,1,ih)
        tta(i,j,1,ih,2)=tt(i,j,1,ih)
        rsa(i,j,1,ih,2)=rs(i,j,1,ih)
      END DO
    END DO

    DO k= 2, ict-1
      DO j= 1, n
        DO i= 1, m
          denm = ts(i,j,k,ih)/( 1.-rsa(i,j,k-1,ih,1)*rs(i,j,k,ih))
          tda(i,j,k,ih,1)= tda(i,j,k-1,ih,1)*td(i,j,k,ih)
          tta(i,j,k,ih,1)= tda(i,j,k-1,ih,1)*tt(i,j,k,ih)               &
                        +(tda(i,j,k-1,ih,1)*rr(i,j,k,ih)                &
                        *rsa(i,j,k-1,ih,1)+tta(i,j,k-1,ih,1))*denm
          rsa(i,j,k,ih,1)= rs(i,j,k,ih)+ts(i,j,k,ih)                    &
                        *rsa(i,j,k-1,ih,1)*denm
          tda(i,j,k,ih,2)= tda(i,j,k,ih,1)
          tta(i,j,k,ih,2)= tta(i,j,k,ih,1)
          rsa(i,j,k,ih,2)= rsa(i,j,k,ih,1)
        END DO
      END DO
    END DO

!-----for middle clouds

    DO im=1,2

      DO k= ict, icb-1
        DO j= 1, n
          DO i= 1, m
            denm = ts(i,j,k,im)/( 1.-rsa(i,j,k-1,ih,im)*rs(i,j,k,im))
            tda(i,j,k,ih,im)= tda(i,j,k-1,ih,im)*td(i,j,k,im)
            tta(i,j,k,ih,im)= tda(i,j,k-1,ih,im)*tt(i,j,k,im)           &
                          +(tda(i,j,k-1,ih,im)*rr(i,j,k,im)             &
                          *rsa(i,j,k-1,ih,im)+tta(i,j,k-1,ih,im))*denm
            rsa(i,j,k,ih,im)= rs(i,j,k,im)+ts(i,j,k,im)                 &
                          *rsa(i,j,k-1,ih,im)*denm
          END DO
        END DO
      END DO

    END DO
  END DO

!-----layers are added one at a time, going up from the surface.
!  rra is the composite reflectance illuminated by beam radiation
!  rxa is the composite reflectance illuminated from above
!      by diffuse radiation
!  rra and rxa are computed from eqs. (4a) and (3a) of Chou

!-----for the low clouds

  DO is=1,2

    DO j= 1, n
      DO i= 1, m
        rra(i,j,np+1,1,is)=rr(i,j,np+1,is)
        rxa(i,j,np+1,1,is)=rs(i,j,np+1,is)
        rra(i,j,np+1,2,is)=rr(i,j,np+1,is)
        rxa(i,j,np+1,2,is)=rs(i,j,np+1,is)
      END DO
    END DO

    DO k=np,icb,-1
      DO j= 1, n
        DO i= 1, m
          denm=ts(i,j,k,is)/( 1.-rs(i,j,k,is)*rxa(i,j,k+1,1,is) )
          rra(i,j,k,1,is)=rr(i,j,k,is)+(td(i,j,k,is)                    &
              *rra(i,j,k+1,1,is)+tt(i,j,k,is)*rxa(i,j,k+1,1,is))*denm
          rxa(i,j,k,1,is)= rs(i,j,k,is)+ts(i,j,k,is)                    &
              *rxa(i,j,k+1,1,is)*denm
          rra(i,j,k,2,is)=rra(i,j,k,1,is)
          rxa(i,j,k,2,is)=rxa(i,j,k,1,is)
        END DO
      END DO
    END DO

!-----for middle clouds

    DO im=1,2

      DO k= icb-1,ict,-1
        DO j= 1, n
          DO i= 1, m
            denm=ts(i,j,k,im)/( 1.-rs(i,j,k,im)*rxa(i,j,k+1,im,is) )
            rra(i,j,k,im,is)= rr(i,j,k,im)+(td(i,j,k,im)                &
                *rra(i,j,k+1,im,is)+tt(i,j,k,im)*rxa(i,j,k+1,im,is))*denm
            rxa(i,j,k,im,is)= rs(i,j,k,im)+ts(i,j,k,im)                 &
                *rxa(i,j,k+1,im,is)*denm
          END DO
        END DO
      END DO

    END DO
  END DO

!-----integration over eight sky situations.
!  ih, im, is denotes high, middle and low cloud groups.

  DO ih=1,2

!-----clear portion

    IF(ih == 1) THEN
      DO j=1,n
        DO i=1,m
          ch(i,j)=1.0-cc(i,j,1)
        END DO
      END DO

    ELSE

!-----cloudy portion

      DO j=1,n
        DO i=1,m
          ch(i,j)=cc(i,j,1)
        END DO
      END DO

    END IF

    DO im=1,2

!-----clear portion

      IF(im == 1) THEN

        DO j=1,n
          DO i=1,m
            cm(i,j)=ch(i,j)*(1.0-cc(i,j,2))
          END DO
        END DO

      ELSE

!-----cloudy portion

        DO j=1,n
          DO i=1,m
            cm(i,j)=ch(i,j)*cc(i,j,2)
          END DO
        END DO

      END IF

      DO is=1,2

!-----clear portion

        IF(is == 1) THEN

          DO j=1,n
            DO i=1,m
              ct(i,j)=cm(i,j)*(1.0-cc(i,j,3))
            END DO
          END DO

        ELSE

!-----cloudy portion

          DO j=1,n
            DO i=1,m
              ct(i,j)=cm(i,j)*cc(i,j,3)
            END DO
          END DO

        END IF

!-----add one layer at a time, going down.

        DO k= icb, np
          DO j= 1, n
            DO i= 1, m
              denm = ts(i,j,k,is)/( 1.-rsa(i,j,k-1,ih,im)*rs(i,j,k,is) )
              tda(i,j,k,ih,im)= tda(i,j,k-1,ih,im)*td(i,j,k,is)
              tta(i,j,k,ih,im)=  tda(i,j,k-1,ih,im)*tt(i,j,k,is)        &
                   +(tda(i,j,k-1,ih,im)*rr(i,j,k,is)                    &
                   *rsa(i,j,k-1,ih,im)+tta(i,j,k-1,ih,im))*denm
              rsa(i,j,k,ih,im)= rs(i,j,k,is)+ts(i,j,k,is)               &
                   *rsa(i,j,k-1,ih,im)*denm
            END DO
          END DO
        END DO

!-----add one layer at a time, going up.

        DO k= ict-1,1,-1
          DO j= 1, n
            DO i= 1, m
              denm =ts(i,j,k,ih)/(1.-rs(i,j,k,ih)*rxa(i,j,k+1,im,is))
              rra(i,j,k,im,is)= rr(i,j,k,ih)+(td(i,j,k,ih)              &
                  *rra(i,j,k+1,im,is)+tt(i,j,k,ih)*rxa(i,j,k+1,im,is))*denm
              rxa(i,j,k,im,is)= rs(i,j,k,ih)+ts(i,j,k,ih)               &
                  *rxa(i,j,k+1,im,is)*denm
            END DO
          END DO
        END DO

!-----compute fluxes following eq (5) of Chou (1992)

!  fdndir is the direct  downward flux
!  fdndif is the diffuse downward flux
!  fupdif is the diffuse upward flux

        DO k=2,np+1
          DO j=1, n
            DO i=1, m
              denm= 1./(1.- rxa(i,j,k,im,is)*rsa(i,j,k-1,ih,im))
              fdndir(i,j)= tda(i,j,k-1,ih,im)
              xx = tda(i,j,k-1,ih,im)*rra(i,j,k,im,is)
              fdndif(i,j)= (xx*rsa(i,j,k-1,ih,im)+tta(i,j,k-1,ih,im))*denm
              fupdif= (xx+tta(i,j,k-1,ih,im)*rxa(i,j,k,im,is))*denm
              flxdn(i,j,k)=fdndir(i,j)+fdndif(i,j)-fupdif
            END DO
          END DO
        END DO

        DO j=1, n
          DO i=1, m
            flxdn(i,j,1)=1.0-rra(i,j,1,im,is)
          END DO
        END DO

!-----summation of fluxes over all (eight) sky situations.

        DO k=1,np+1
          DO j=1,n
            DO i=1,m
              IF(ih == 1 .AND. im == 1 .AND. is == 1) THEN
                fclr(i,j,k)=flxdn(i,j,k)
              END IF
              fall(i,j,k)=fall(i,j,k)+flxdn(i,j,k)*ct(i,j)
            END DO
          END DO
        END DO

        DO j=1,n
          DO i=1,m
            fsdir(i,j)=fsdir(i,j)+fdndir(i,j)*ct(i,j)
            fsdif(i,j)=fsdif(i,j)+fdndif(i,j)*ct(i,j)
          END DO
        END DO

      END DO
    END DO
  END DO

  RETURN
END SUBROUTINE cldflx