pkg/fizhi/fizhi_swrad.F

C $Header: /u/gcmpack/MITgcm/pkg/fizhi/fizhi_swrad.F,v 1.24 2005/06/16 16:46:12 ce107 Exp $
C $Name:  $

#include "FIZHI_OPTIONS.h"
      subroutine swrio (nymd,nhms,bi,bj,ndswr,myid,istrip,npcs,
     .        low_level,mid_level,im,jm,lm,
     .        pz,plz,plze,dpres,pkht,pkz,tz,qz,oz,co2,
     .        albvisdr,albvisdf,albnirdr,albnirdf,
     .        dtradsw,dtswclr,radswg,swgclr,
     .        fdifpar,fdirpar,osr,osrclr,
     .        ptop,nswcld,cldsw,cswmo,nswlz,swlz,
     .        lpnt,imstturb,qliqave,fccave,landtype,xlats,xlons)

      implicit none

c Input Variables
c ---------------
      integer nymd,nhms,bi,bj,ndswr,myid,istrip,npcs
      integer mid_level,low_level
      integer im,jm,lm        
      _RL  ptop
      _RL pz(im,jm),plz(im,jm,lm),plze(im,jm,lm+1),dpres(im,jm,lm)
      _RL pkht(im,jm,lm+1),pkz(im,jm,lm)
      _RL tz(im,jm,lm),qz(im,jm,lm)
      _RL oz(im,jm,lm)
      _RL co2
      _RL albvisdr(im,jm),albvisdf(im,jm),albnirdr(im,jm)
      _RL albnirdf(im,jm)
      _RL radswg(im,jm),swgclr(im,jm),fdifpar(im,jm),fdirpar(im,jm)
      _RL osr(im,jm),osrclr(im,jm),dtradsw(im,jm,lm)
      _RL dtswclr(im,jm,lm)
      integer nswcld,nswlz    
      _RL cldsw(im,jm,lm),cswmo(im,jm,lm),swlz(im,jm,lm)   
      logical lpnt            
      integer imstturb        
      _RL qliqave(im,jm,lm),fccave(im,jm,lm)  
      integer landtype(im,jm) 
      _RL xlats(im,jm),xlons(im,jm)

c Local Variables
c ---------------
      integer   i,j,L,nn,nsecf
      integer   ntmstp,nymd2,nhms2
      _RL      getcon,grav,cp,undef
      _RL      ra,alf,reffw,reffi,tminv

      parameter ( reffw = 10.0 )   
      parameter ( reffi = 65.0 )   

      _RL tdry(im,jm,lm)
      _RL TEMP1(im,jm)
      _RL TEMP2(im,jm)
      _RL zenith (im,jm)
      _RL cldtot (im,jm,lm)
      _RL cldmxo (im,jm,lm)
      _RL totcld (im,jm)
      _RL cldlow (im,jm)
      _RL cldmid (im,jm)
      _RL cldhi  (im,jm)
      _RL taulow (im,jm)
      _RL taumid (im,jm)
      _RL tauhi  (im,jm)
      _RL tautype(im,jm,lm,3)
      _RL tau(im,jm,lm)
      _RL albedo(im,jm)    

      _RL PK(ISTRIP,lm)
      _RL qzl(ISTRIP,lm),CLRO(ISTRIP,lm)
      _RL TZL(ISTRIP,lm)
      _RL OZL(ISTRIP,lm)
      _RL PLE(ISTRIP,lm+1)
      _RL COSZ(ISTRIP)
      _RL dpstrip(ISTRIP,lm)

      _RL albuvdr(ISTRIP),albuvdf(ISTRIP)
      _RL albirdr(ISTRIP),albirdf(ISTRIP)
      _RL difpar (ISTRIP),dirpar (ISTRIP)

      _RL fdirir(istrip),fdifir(istrip)
      _RL fdiruv(istrip),fdifuv(istrip)

      _RL flux(istrip,lm+1)
      _RL fluxclr(istrip,lm+1)
      _RL dtsw(istrip,lm)
      _RL dtswc(istrip,lm)

      _RL taul   (istrip,lm)
      _RL reff   (istrip,lm,2)
      _RL tauc   (istrip,lm,2)
      _RL taua   (istrip,lm)
      _RL tstrip (istrip)
#ifdef ALLOW_DIAGNOSTICS
      logical  diagnostics_is_on
      external diagnostics_is_on
      _RL tmpdiag(im,jm),tmpdiag2(im,jm)
#endif

      logical first
      data first /.true./

C **********************************************************************
C ****                       INITIALIZATION                         ****
C **********************************************************************

      grav  = getcon('GRAVITY')
      cp    = getcon('CP')
      undef = getcon('UNDEF')

      NTMSTP = nsecf(NDSWR)
      TMINV  = 1./float(ntmstp)

C Compute Temperature from Theta
C ------------------------------
      do L=1,lm
      do j=1,jm
      do i=1,im
      tdry(i,j,L) = tz(i,j,L)*pkz(i,j,L)
      enddo
      enddo
      enddo

      if (first .and. myid.eq.1 ) then
      print *
      print *,'Low-Level Clouds are Grouped between levels: ',
     .         lm,' and ',low_level
      print *,'Mid-Level Clouds are Grouped between levels: ',
     .         low_level-1,' and ',mid_level
      print *
      first = .false.
      endif

C **********************************************************************
C ****             CALCULATE COSINE OF THE ZENITH ANGLE             ****
C **********************************************************************

      CALL ASTRO ( NYMD,   NHMS,  XLATS,XLONS, im*jm, TEMP1,RA )
                   NYMD2 = NYMD
                   NHMS2 = NHMS
      CALL TICK  ( NYMD2,  NHMS2, NTMSTP )
      CALL ASTRO ( NYMD2,  NHMS2, XLATS,XLONS, im*jm, TEMP2,RA )

      do j = 1,jm
      do i = 1,im
          zenith(I,j) = TEMP1(I,j) + TEMP2(I,j)
      IF( zenith(I,j) .GT. 0.0 ) THEN
          zenith(I,j) = (2./3.)*( zenith(I,j)-TEMP2(I,j)*TEMP1(I,j)
     .                          / zenith(I,j) )
      ENDIF
      ENDDO
      ENDDO

C **********************************************************************
c ****        Compute Two-Dimension Total Cloud Fraction (0-1)      ****
C **********************************************************************

c Initialize Clear Sky Probability for Low, Mid, and High Clouds
c --------------------------------------------------------------
      do j =1,jm
      do i =1,im
      cldlow(i,j) = 0.0
      cldmid(i,j) = 0.0
      cldhi (i,j) = 0.0
      enddo
      enddo

c Adjust cloud fractions and cloud liquid water due to moist turbulence
c ---------------------------------------------------------------------
      if(imstturb.ne.0) then
        do L =1,lm
        do j =1,jm
        do i =1,im
           cldtot(i,j,L)=min(1.0 _d 0,max(cldsw(i,j,L),fccave(i,j,L)/
     $          imstturb))
           cldmxo(i,j,L)=min(1.0 _d 0,cswmo(i,j,L))
           swlz(i,j,L)=swlz(i,j,L)+qliqave(i,j,L)/imstturb
        enddo
        enddo
        enddo
      else
        do L =1,lm
        do j =1,jm
        do i =1,im
         cldtot(i,j,L) =  min( 1.0 _d 0,cldsw(i,j,L) )
         cldmxo(i,j,L) =  min( 1.0 _d 0,cswmo(i,j,L) )
        enddo
        enddo
        enddo
      endif

c Compute 3-D Cloud Fractions
c ---------------------------
      if( nswcld.ne.0 ) then
      do L = 1,lm
      do j = 1,jm
      do i = 1,im
c Compute Low-Mid-High Maximum Overlap Cloud Fractions
c ----------------------------------------------------
      if(     L.lt.mid_level ) then
             cldhi (i,j) = max( cldtot(i,j,L),cldhi (i,j) )
      elseif( L.ge.mid_level .and. 
     .        L.lt.low_level ) then
             cldmid(i,j) = max( cldtot(i,j,L),cldmid(i,j) )
      elseif( L.ge.low_level ) then
             cldlow(i,j) = max( cldtot(i,j,L),cldlow(i,j) )
      endif

      enddo
      enddo
      enddo
      endif

c Totcld => Product of Clear Sky Probabilities for Low, Mid, and High Clouds
c --------------------------------------------------------------------------
      do j = 1,jm
      do i = 1,im
      totcld(i,j) = 1.0 - (1.-cldhi (i,j))
     .                  * (1.-cldmid(i,j))
     .                  * (1.-cldlow(i,j))
      enddo
      enddo

#ifdef ALLOW_DIAGNOSTICS

c Compute Cloud Diagnostics
c -------------------------
      if(diagnostics_is_on('CLDFRC  ',myid) ) then
       call diagnostics_fill(totcld,'CLDFRC  ',0,1,3,bi,bj,myid)
      endif

      if(diagnostics_is_on('CLDRAS  ',myid) ) then
       do L=1,lm
       do j=1,jm
       do i=1,im
        tmpdiag(i,j) = cswmo(i,j,L)
       enddo
       enddo
       call diagnostics_fill(tmpdiag,'CLDRAS  ',L,1,3,bi,bj,myid)
      enddo
      endif

      if(diagnostics_is_on('CLDTOT  ',myid) ) then
       do L=1,lm
       do j=1,jm
       do i=1,im
        tmpdiag(i,j) = cldtot(i,j,L)
       enddo
       enddo
       call diagnostics_fill(tmpdiag,'CLDTOT  ',L,1,3,bi,bj,myid)
      enddo
      endif

      if(diagnostics_is_on('CLDLOW  ',myid) ) then
       call diagnostics_fill(cldlow,'CLDLOW  ',0,1,3,bi,bj,myid)
      endif

      if(diagnostics_is_on('CLDMID  ',myid) ) then
       call diagnostics_fill(cldmid,'CLDMID  ',0,1,3,bi,bj,myid)
      endif

      if(diagnostics_is_on('CLDHI   ',myid) ) then
       call diagnostics_fill(cldhi,'CLDHI   ',0,1,3,bi,bj,myid)
      endif

      if(diagnostics_is_on('LZRAD   ',myid) ) then
       do L=1,lm
       do j=1,jm
       do i=1,im
        tmpdiag(i,j) = swlz(i,j,L) * 1.0e6
       enddo
       enddo
       call diagnostics_fill(tmpdiag,'LZRAD   ',L,1,3,bi,bj,myid)
       enddo
      endif

c Albedo Diagnostics
c ------------------
      if(diagnostics_is_on('ALBVISDR',myid) ) then
       call diagnostics_fill(albvisdr,'ALBVISDR',0,1,3,bi,bj,myid)
      endif

      if(diagnostics_is_on('ALBVISDF',myid) ) then
       call diagnostics_fill(albvisdf,'ALBVISDF',0,1,3,bi,bj,myid)
      endif

      if(diagnostics_is_on('ALBNIRDR',myid) ) then
       call diagnostics_fill(albnirdr,'ALBNIRDR',0,1,3,bi,bj,myid)
      endif

      if(diagnostics_is_on('ALBNIRDF',myid) ) then
       call diagnostics_fill(albnirdf,'ALBNIRDF',0,1,3,bi,bj,myid)
      endif

#endif

C Compute Optical Thicknesses and Diagnostics
C -------------------------------------------
      call opthk(tdry,plz,plze,swlz,cldtot,cldmxo,landtype,im,jm,lm,
     .                                                          tautype)

      do L = 1,lm
      do j = 1,jm
      do i = 1,im
      tau(i,j,L) = tautype(i,j,L,1)+tautype(i,j,L,2)+tautype(i,j,L,3)
      enddo
      enddo
      enddo

#ifdef ALLOW_DIAGNOSTICS
      if(diagnostics_is_on('TAUAVE  ',myid) ) then
       do L=1,lm
       do j=1,jm
       do i=1,im
        tmpdiag(i,j) = tau(i,j,L)*100/(plze(i,j,L+1)-plze(i,j,L))
       enddo
       enddo
       call diagnostics_fill(tmpdiag,'TAUAVE  ',L,1,3,bi,bj,myid)
      enddo
      endif

      if(diagnostics_is_on('TAUCLD  ',myid) .and.
     .                 diagnostics_is_on('TAUCLDC ',myid) ) then
       do L=1,lm
       do j=1,jm
       do i=1,im
        if( cldtot(i,j,L).ne.0.0 ) then
        tmpdiag(i,j) = tau(i,j,L)*100/(plze(i,j,L+1)-plze(i,j,L))
        tmpdiag2(i,j) = 1.
        else
        tmpdiag(i,j) = 0.
        tmpdiag2(i,j) = 0.
        endif
       enddo
       enddo
       call diagnostics_fill(tmpdiag,'TAUCLD  ',L,1,3,bi,bj,myid)
       call diagnostics_fill(tmpdiag,'TAUCLDC ',L,1,3,bi,bj,myid)
      enddo
      endif

c Compute Low, Mid, and High Cloud Optical Depth Diagnostics
c ----------------------------------------------------------
      if(diagnostics_is_on('TAULOW  ',myid) .and.
     .                 diagnostics_is_on('TAULOWC ',myid) ) then
       do j = 1,jm
       do i = 1,im
        taulow(i,j) =  0.0
        if( cldlow(i,j).ne.0.0 ) then
         do L = low_level,lm
          taulow(i,j) = taulow(i,j) + tau(i,j,L)
         enddo
         tmpdiag2(i,j) = 1.
        else
         tmpdiag(i,j) = 0.
        endif
       enddo
       enddo
       call diagnostics_fill(taulow,'TAULOW  ',0,1,3,bi,bj,myid)
       call diagnostics_fill(tmpdiag2,'TAULOWC ',0,1,3,bi,bj,myid)
      endif

      if(diagnostics_is_on('TAUMID  ',myid) .and.
     .                 diagnostics_is_on('TAUMIDC ',myid) ) then
       do j = 1,jm
       do i = 1,im
        taumid(i,j) =  0.0
        if( cldmid(i,j).ne.0.0 ) then
         do L = mid_level,low_level+1
          taumid(i,j) = taumid(i,j) + tau(i,j,L)
         enddo
         tmpdiag2(i,j) = 1.
        else
         tmpdiag(i,j) = 0.
        endif
       enddo
       enddo
       call diagnostics_fill(taumid,'TAUMID  ',0,1,3,bi,bj,myid)
       call diagnostics_fill(tmpdiag2,'TAUMIDC ',0,1,3,bi,bj,myid)
      endif

      if(diagnostics_is_on('TAUHI   ',myid) .and.
     .                 diagnostics_is_on('TAUHIC  ',myid) ) then
       do j = 1,jm
       do i = 1,im
        tauhi(i,j) =  0.0
        if( cldhi(i,j).ne.0.0 ) then
         do L = 1,mid_level+1
          tauhi(i,j) = tauhi(i,j) + tau(i,j,L)
         enddo
         tmpdiag2(i,j) = 1.
        else
         tmpdiag(i,j) = 0.
        endif
       enddo
       enddo
       call diagnostics_fill(tauhi,'TAUHI   ',0,1,3,bi,bj,myid)
       call diagnostics_fill(tmpdiag2,'TAUHIC  ',0,1,3,bi,bj,myid)
      endif
#endif

C***********************************************************************
C ****                     LOOP OVER GLOBAL REGIONS                 ****
C **********************************************************************

      do 1000 nn = 1,npcs

C **********************************************************************
C ****                   VARIABLE INITIALIZATION                    ****
C **********************************************************************

      CALL STRIP ( zenith,COSZ,im*jm,ISTRIP,1,NN )

      CALL STRIP ( plze,  ple   ,im*jm,ISTRIP,lm+1,NN)
      CALL STRIP ( pkz ,  pk    ,im*jm,ISTRIP,lm  ,NN)
      CALL STRIP ( dpres,dpstrip,im*jm,ISTRIP,lm  ,NN)
      CALL STRIP ( tdry,  tzl   ,im*jm,ISTRIP,lm  ,NN)
      CALL STRIP ( qz  ,  qzl   ,im*jm,ISTRIP,lm  ,NN)
      CALL STRIP ( oz  ,  ozl   ,im*jm,ISTRIP,lm  ,NN)
      CALL STRIP ( tau ,  taul  ,im*jm,ISTRIP,lm  ,NN)

      CALL STRIP ( albvisdr,albuvdr,im*jm,ISTRIP,1,NN )
      CALL STRIP ( albvisdf,albuvdf,im*jm,ISTRIP,1,NN )
      CALL STRIP ( albnirdr,albirdr,im*jm,ISTRIP,1,NN )
      CALL STRIP ( albnirdf,albirdf,im*jm,ISTRIP,1,NN )

      call strip ( cldtot,clro,im*jm,istrip,lm,nn )

c Partition Tau between Water and Ice Particles
c ---------------------------------------------
      do L= 1,lm
      do i= 1,istrip
              alf = min( max((tzl(i,l)-253.15)/20.,0. _d 0) ,1. _d 0)
      taua(i,L)   = 0.

      if( alf.ne.0.0 .and. alf.ne.1.0 ) then
      tauc(i,L,1) = taul(i,L)/(1.+alf/(1-alf) * (reffi/reffw*0.80) )
      tauc(i,L,2) = taul(i,L)/(1.+(1-alf)/alf * (reffw/reffi*1.25) )
      endif

      if( alf.eq.0.0 ) then
      tauc(i,L,1) = taul(i,L)
      tauc(i,L,2) = 0.0
      endif

      if( alf.eq.1.0 ) then
      tauc(i,L,1) = 0.0
      tauc(i,L,2) = taul(i,L)
      endif

      reff(i,L,1) = reffi
      reff(i,L,2) = reffw
      enddo
      enddo

      call sorad ( istrip,1,1,lm,ple,tzl,qzl,ozl,co2,
     .             tauc,reff,clro,mid_level,low_level,taua,
     .             albirdr,albirdf,albuvdr,albuvdf,cosz,
     .             flux,fluxclr,fdirir,fdifir,dirpar,difpar,
     .             fdiruv,fdifuv )

C **********************************************************************
C ****     Compute Mass-Weighted Theta Tendencies from Fluxes       ****
C **********************************************************************

      do l=1,lm
      do i=1,istrip
      alf = grav*(ple(i,Lm+1)-ptop)/(cp*dpstrip(i,L)*100)
      dtsw (i,L) = alf*( flux   (i,L)-flux   (i,L+1) )/pk(i,L)
      dtswc(i,L) = alf*( fluxclr(i,L)-fluxclr(i,L+1) )/pk(i,L)
      enddo
      enddo
      
      call paste ( dtsw , dtradsw ,istrip,im*jm,lm,nn )
      call paste ( dtswc, dtswclr ,istrip,im*jm,lm,nn )

      call paste ( flux   (1,1),osr   ,istrip,im*jm,1,nn )
      call paste ( fluxclr(1,1),osrclr,istrip,im*jm,1,nn )

      call paste ( flux   (1,lm+1),radswg,istrip,im*jm,1,nn )
      call paste ( fluxclr(1,lm+1),swgclr,istrip,im*jm,1,nn )

      call paste ( dirpar,fdirpar,istrip,im*jm,1,nn )
      call paste ( difpar,fdifpar,istrip,im*jm,1,nn )

c Calculate Mean Albedo
c ---------------------
      do i=1,istrip
       if( cosz(i).gt.0.0 ) then
        tstrip(i) = 1.0 - flux(i,lm+1)/
     . ( fdirir(i)+fdifir(i)+dirpar(i)+difpar(i) + fdiruv(i)+fdifuv(i) )
        if( tstrip(i).lt.0.0 ) tstrip(i) = undef
       else
        tstrip(i) = undef
       endif
      enddo
      call paste ( tstrip,albedo,istrip,im*jm,1,nn )

 1000 CONTINUE

#ifdef ALLOW_DIAGNOSTICS

c Mean Albedo Diagnostic
c ----------------------
      if(diagnostics_is_on('ALBEDO  ',myid) .and.
     .                 diagnostics_is_on('ALBEDOC ',myid) ) then
       do j=1,jm
       do i=1,im
        if( albedo(i,j).ne.undef ) then
         tmpdiag(i,j) = albedo(i,j)
         tmpdiag2(i,j) = 1.
        else
         tmpdiag(i,j) = 0.
         tmpdiag2(i,j) = 0.
        endif
       enddo
       enddo
       call diagnostics_fill(tmpdiag,'ALBEDO  ',0,1,3,bi,bj,myid)
       call diagnostics_fill(tmpdiag2,'ALBEDOC ',0,1,3,bi,bj,myid)
      endif
#endif

C **********************************************************************
C ****                 ZERO-OUT OR BUMP COUNTERS                    ****
C **********************************************************************

      nswlz    = 0
      nswcld   = 0
      imstturb = 0

      do L = 1,lm
      do j = 1,jm
      do i = 1,im
         fccave(i,j,L) = 0.
        qliqave(i,j,L) = 0.
      enddo
      enddo
      enddo

      return
      end
      subroutine opthk ( tl,pl,ple,lz,cf,cfm,lwi,im,jm,lm,tau )
C***********************************************************************
C
C  PURPOSE:
C  ========
C    Compute cloud optical thickness using temperature and
C    cloud fractions.
C
C  INPUT:
C  ======
C    tl ......... Temperature at Model Layers (K)
C    pl ......... Model Layer Pressures (mb)
C    ple ........ Model Edge  Pressures (mb)
C    lz ......... Diagnosed Convective and Large-Scale Cloud Water (g/g)
C    cf ......... Total Cloud Fraction (Random + Maximum Overlap)
C    cfm ........ Maximum Overlap Cloud Fraction
C    lwi ........ Surface Land Type
C    im ......... Longitudinal dimension
C    jm ......... Latitudinal  dimension
C    lm ......... Vertical     dimension
C
C  OUTPUT:
C  =======
C    tau ........ Cloud Optical Thickness (non-dimensional)
C                  tau(im,jm,lm,1):  Suspended Ice
C                  tau(im,jm,lm,2):  Suspended Water
C                  tau(im,jm,lm,3):  Raindrops
C
C***********************************************************************

      implicit none

      integer  im,jm,lm,i,j,L

      _RL  tl(im,jm,lm)
      _RL  pl(im,jm,lm)
      _RL ple(im,jm,lm+1)
      _RL  lz(im,jm,lm)
      _RL  cf(im,jm,lm)
      _RL cfm(im,jm,lm)
      _RL tau(im,jm,lm,3)
      integer lwi(im,jm)

      _RL dp, alf, fracls, fraccu
      _RL tauice, tauh2o, tauras

c Compute Cloud Optical Depths
c ----------------------------
      do L=1,lm
      do j=1,jm
      do i=1,im
         alf   =  min( max(( tl(i,j,L)-233.15)/20.,0. _d 0), 1. _d 0)
         dp   =  ple(i,j,L+1)-ple(i,j,L)
         tau(i,j,L,1)  = 0.0
         tau(i,j,L,2)  = 0.0
         tau(i,j,L,3)  = 0.0
      if( cf(i,j,L).ne.0.0 ) then

c Determine fraction of large-scale and cumulus clouds
c ----------------------------------------------------
         fracls = ( cf(i,j,L)-cfm(i,j,L) )/cf(i,j,L)
         fraccu = 1.0-fracls

c Define tau for large-scale ice and water clouds
c Note:  tauice is scaled between (0.02 & .2) for:  0 < lz <  2 mg/kg over Land 
c Note:  tauh2o is scaled between (0.20 & 20) for:  0 < lz <  5 mg/kg over Land 
c Note:  tauh2o is scaled between (0.20 & 12) for:  0 < lz < 50 mg/kg over Ocean
c -------------------------------------------------------------------------------

c Large-Scale Ice
c ---------------
         tauice = max( 0.0002 _d 0, 0.002*min(500*lz(i,j,L)*1000,
     $        1.0 _d 0) )
         tau(i,j,L,1) = fracls*(1-alf)*tauice*dp

c Large-Scale Water
c -----------------
C Over Land
         if( lwi(i,j).le.10 ) then
            tauh2o = max( 0.0020 _d 0, 0.200*min(200*lz(i,j,L)*1000,
     $           1.0 _d 0) )  
          tau(i,j,L,3) = fracls*alf*tauh2o*dp
         else
C Non-Precipitation Clouds Over Ocean
          if( lz(i,j,L).eq.0.0 ) then
           tauh2o = .12                      
           tau(i,j,L,2) = fracls*alf*tauh2o*dp
          else
C Over Ocean
             tauh2o = max( 0.0020 _d 0, 0.120*min( 20*lz(i,j,L)*1000,
     $            1.0 _d 0) )
           tau(i,j,L,3) = fracls*alf*tauh2o*dp
          endif
         endif

c Sub-Grid Convective
c -------------------
         if( tl(i,j,L).gt.210.0 ) then
         tauras = .16
         tau(i,j,L,3) = tau(i,j,L,3) + fraccu*tauras*dp
         else
         tauras = tauice
         tau(i,j,L,1) = tau(i,j,L,1) + fraccu*tauras*dp
         endif

c Define tau for large-scale and cumulus clouds
c ---------------------------------------------
ccc      tau(i,j,L) = ( fracls*( alf*tauh2o + (1-alf)*tauice )
ccc  .                + fraccu*tauras )*dp

      endif

      enddo
      enddo
      enddo

      return
      end
      subroutine sorad(m,n,ndim,np,pl,ta,wa,oa,co2,
     *                 taucld,reff,fcld,ict,icb,
     *                 taual,rsirbm,rsirdf,rsuvbm,rsuvdf,cosz,
     *                 flx,flc,fdirir,fdifir,fdirpar,fdifpar,
     *                 fdiruv,fdifuv)
c********************************************************************
c  
c This routine computes solar fluxes due to the absoption by water
c  vapor, ozone, co2, o2, clouds, and aerosols and due to the
c  scattering by clouds, aerosols, and gases.
c
c This is a vectorized code. It computes the fluxes simultaneous for
c  (m x n) soundings, which is a subset of the (m x ndim) soundings.
c  In a global climate model, m and ndim correspond to the numbers of
c  grid boxes in the zonal and meridional directions, respectively.
c
c Ice and liquid cloud particles are allowed to co-exist in any of the
c  np layers. Two sets of cloud parameters are required as inputs, one
c  for ice paticles and the other for liquid particles.  These parameters
c  are optical thickness (taucld) and effective particle size (reff).
c
c If no information is available for reff, a default value of
c  10 micron for liquid water and 75 micron for ice can be used.
c
c Clouds are grouped into high, middle, and low clouds separated by the
c  level indices ict and icb.  For detail, see the subroutine cldscale.
c
c----- Input parameters:                           
c                                                   units      size
c   number of soundings in zonal direction (m)       n/d        1
c   number of soundings in meridional direction (n)  n/d        1
c   maximum number of soundings in                   n/d        1
c           meridional direction (ndim)
c   number of atmospheric layers (np)                n/d        1
c   level pressure (pl)                              mb       m*ndim*(np+1)
c   layer temperature (ta)                           k        m*ndim*np
c   layer specific humidity (wa)                     gm/gm    m*ndim*np
c   layer ozone concentration (oa)                   gm/gm    m*ndim*np
c   co2 mixing ratio by volumn (co2)               parts/part   1
c   cloud optical thickness (taucld)                 n/d      m*ndim*np*2
c                index 1 for ice particles
c                index 2 for liquid drops
c   effective cloud-particle size (reff)           micrometer m*ndim*np*2
c                index 1 for ice particles
c                index 2 for liquid drops
c   cloud amount (fcld)                            fraction   m*ndim*np
c   level index separating high and middle           n/d        1
c                clouds (ict)
c   level index separating middle and low clouds     n/d        1
c                clouds (icb)
c   aerosol optical thickness (taual)                n/d      m*ndim*np 
c   solar ir surface albedo for beam                fraction   m*ndim
c                radiation (rsirbm)                
c   solar ir surface albedo for diffuse             fraction   m*ndim
c                radiation (rsirdf)                
c   uv + par surface albedo for beam                     fraction   m*ndim
c                radiation (rsuvbm)                
c   uv + par surface albedo for diffuse                  fraction   m*ndim
c                radiation (rsuvdf)                
c   cosine of solar zenith angle (cosz)            n/d        m*ndim
c
c----- Output parameters
c
c   all-sky flux (downward minus upward) (flx)     fraction   m*ndim*(np+1)
c   clear-sky flux (downward minus upward) (flc)   fraction   m*ndim*(np+1)
c   all-sky direct downward ir (0.7-10 micron)
c                flux at the surface (fdirir)      fraction   m*ndim
c   all-sky diffuse downward ir flux at
c                the surface (fdifir)              fraction   m*ndim
c   all-sky direct downward par (0.4-0.7 micron)
c                flux at the surface (fdirpar)     fraction   m*ndim
c   all-sky diffuse downward par flux at
c                the surface (fdifpar)             fraction   m*ndim
*
c----- Notes:
c
c    (1) The unit of flux is fraction of the incoming solar radiation
c        at the top of the atmosphere.  Therefore, fluxes should
c        be equal to flux multiplied by the extra-terrestrial solar
c        flux and the cosine of solar zenith angle.
c    (2) Clouds and aerosols can be included in any layers by specifying
c        fcld(i,j,k), taucld(i,j,k,*) and taual(i,j,k), k=1,np. 
c        For an atmosphere without clouds and aerosols,
c        set fcld(i,j,k)=taucld(i,j,k,*)=taual(i,j,k)=0.0.
c    (3) Aerosol single scattering albedos and asymmetry
c        factors are specified in the subroutines solir and soluv.
c    (4) pl(i,j,1) is the pressure at the top of the model, and
c        pl(i,j,np+1) is the surface pressure.
c    (5) the pressure levels ict and icb correspond approximately
c        to 400 and 700 mb.
c        
c**************************************************************************
 
      implicit none

c-----Explicit Inline Directives

#ifdef CRAY
#ifdef f77
cfpp$ expand (expmn)
#endif
#endif
      _RL expmn

c-----input parameters

      integer m,n,ndim,np,ict,icb
      _RL pl(m,ndim,np+1),ta(m,ndim,np),wa(m,ndim,np),oa(m,ndim,np)
      _RL  taucld(m,ndim,np,2),reff(m,ndim,np,2)
      _RL  fcld(m,ndim,np),taual(m,ndim,np)
      _RL  rsirbm(m,ndim),rsirdf(m,ndim),
     *     rsuvbm(m,ndim),rsuvdf(m,ndim),cosz(m,ndim),co2

c-----output parameters

      _RL  flx(m,ndim,np+1),flc(m,ndim,np+1)
      _RL  fdirir(m,ndim),fdifir(m,ndim)
      _RL  fdirpar(m,ndim),fdifpar(m,ndim)
      _RL  fdiruv(m,ndim),fdifuv(m,ndim)

c-----temporary array
 
      integer i,j,k
      _RL  cc(m,n,3),tauclb(m,n,np),tauclf(m,n,np)
      _RL  dp(m,n,np),wh(m,n,np),oh(m,n,np),scal(m,n,np)
      _RL  swh(m,n,np+1),so2(m,n,np+1),df(m,n,np+1)
      _RL  sdf(m,n),sclr(m,n),csm(m,n),x
 
c-----------------------------------------------------------------

      do j= 1, n 
       do i= 1, m 

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

c-----csm is the effective secant of the solar zenith angle
c     see equation (12) of Lacis and Hansen (1974, JAS)    
 
         csm(i,j)=35./sqrt(1224.*cosz(i,j)*cosz(i,j)+1.)

       enddo 
      enddo


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

c-----compute layer thickness and pressure-scaling function. 
c     indices for the surface level and surface layer
c     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
 
c-----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)

c-----compute ozone amount, unit is (cm-atm)stp.
 
          oh(i,j,k)=1.02*oa(i,j,k)*dp(i,j,k)*466.7

        enddo
       enddo
      enddo


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

      call cldscale(m,n,ndim,np,cosz,fcld,taucld,ict,icb,
     *              cc,tauclb,tauclf)

c-----initialize fluxes for all-sky (flx), clear-sky (flc), and
c     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.
        enddo
       enddo
      enddo

c-----compute solar ir fluxes

      call solir (m,n,ndim,np,wh,taucld,tauclb,tauclf,reff,ict,icb
     *     ,fcld,cc,taual,csm,rsirbm,rsirdf,flx,flc,fdirir,fdifir)

c-----compute and update uv and par fluxes

      call soluv (m,n,ndim,np,oh,dp,taucld,tauclb,tauclf,reff,ict,icb
     *     ,fcld,cc,taual,csm,rsuvbm,rsuvdf,flx,flc
     *     ,fdirpar,fdifpar,fdiruv,fdifuv)

c-----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)
        enddo
       enddo
      enddo

c-----compute flux reduction due to oxygen following
c      chou (J. climate, 1990). The fraction 0.0287 is the
c      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)))
         enddo
        enddo
       enddo          

c-----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)
        enddo
       enddo
      enddo

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

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

c-----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)
        enddo
       enddo
      enddo

c-----adjust for the all-sky fluxes due to o2 and co2.  It is
c     assumed that o2 and co2 have no effects on solar radiation
c     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
       enddo
      enddo

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

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

         if(fcld(i,j,k).gt.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))
         endif
          flx(i,j,k+1)=flx(i,j,k+1)-sdf(i,j)-df(i,j,k+1)*sclr(i,j)

        enddo
       enddo
      enddo

c-----adjust for the direct downward ir flux.
      do j= 1, n
       do i= 1, m
         x = (1.-rsirdf(i,j))*fdifir(i,j) + (1.-rsirbm(i,j))*fdirir(i,j)
         x = (-sdf(i,j)-df(i,j,np+1)*sclr(i,j))/x
         fdirir(i,j)=fdirir(i,j)*(1.+x)
         fdifir(i,j)=fdifir(i,j)*(1.+x)
       enddo
      enddo

      return
      end  

c********************************************************************

      subroutine cldscale (m,n,ndim,np,cosz,fcld,taucld,ict,icb,
     *                     cc,tauclb,tauclf)

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

      implicit none

c-----input parameters

      integer m,n,ndim,np,ict,icb
      _RL  cosz(m,ndim),fcld(m,ndim,np),taucld(m,ndim,np,2)

c-----output parameters

      _RL  cc(m,n,3),tauclb(m,n,np),tauclf(m,n,np)

c-----temporary variables

      integer i,j,k,im,it,ia,kk
      _RL   fm,ft,fa,xai,taux

c-----pre-computed table

      integer   nm,nt,na
      parameter (nm=11,nt=9,na=11) 
      _RL   dm,dt,da,t1,caib(nm,nt,na),caif(nt,na)
      parameter (dm=0.1,dt=0.30103,da=0.1,t1=-0.9031)

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

#include "cai-dat.h"
      save caib,caif


c-----clouds within each of the high, middle, and low clouds are
c     assumed maximally overlapped, and the cloud cover (cc)
c     for a group is the maximum cloud cover of all the layers
c     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
       enddo
      enddo

      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))
        enddo
       enddo
      enddo

      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))
        enddo
       enddo
      enddo

      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))
        enddo
       enddo
      enddo

c-----scale the cloud optical thickness.
c     taucld(i,j,k,1) is the optical thickness for ice particles, and
c     taucld(i,j,k,2) is the optical thickness for liquid particles.
      
      do k=1,np

         if(k.lt.ict) then
            kk=1
         elseif(k.ge.ict .and. k.lt.icb) then
            kk=2
         else
            kk=3
         endif

       do j=1,n
        do i=1,m

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

         taux=taucld(i,j,k,1)+taucld(i,j,k,2)
         if (taux.gt.0.05 .and. fcld(i,j,k).gt.0.01) then

c-----normalize cloud cover

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

c-----table look-up

           taux=min(taux,32. _d 0)

           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)

c-----scale cloud optical thickness for beam radiation.
c     the scaling factor, xai, is a function of the solar zenith
c     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 _d 0)
     
           tauclb(i,j,k) = taux*xai

c-----scale cloud optical thickness for diffuse radiation.
c     the scaling factor, xai, is a function of the cloud optical
c     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 _d 0)
     
           tauclf(i,j,k) = taux*xai

         endif

        enddo
       enddo
      enddo

      return
      end
c***********************************************************************

      subroutine solir (m,n,ndim,np,wh,taucld,tauclb,tauclf,reff,
     *  ict,icb,fcld,cc,taual,csm,rsirbm,rsirdf,flx,flc,fdirir,fdifir)
 
c************************************************************************
c  compute solar flux in the infrared region. The spectrum is divided
c   into three bands:
c
c          band   wavenumber(/cm)  wavelength (micron)
c           1       1000-4400         2.27-10.0
c           2       4400-8200         1.22-2.27
c           3       8200-14300        0.70-1.22
c
c----- Input parameters:                            units      size
c
c   number of soundings in zonal direction (m)       n/d        1
c   number of soundings in meridional direction (n)  n/d        1
c   maximum number of soundings in                   n/d        1
c          meridional direction (ndim)
c   number of atmospheric layers (np)                n/d        1
c   layer water vapor content (wh)                 gm/cm^2    m*n*np
c   cloud optical thickness (taucld)                 n/d      m*ndim*np*2
c          index 1 for ice paticles
c          index 2 for liquid particles
c   scaled cloud optical thickness                   n/d      m*n*np
c          for beam radiation (tauclb)
c   scaled cloud optical thickness                   n/d      m*n*np
c          for diffuse radiation  (tauclf)
c   effective cloud-particle size (reff)           micrometer m*ndim*np*2
c          index 1 for ice paticles
c          index 2 for liquid particles
c   level index separating high and                  n/d      m*n
c          middle clouds (ict)
c   level index separating middle and                n/d      m*n
c          low clouds (icb)
c   cloud amount (fcld)                            fraction   m*ndim*np
c   cloud amount of high, middle, and                n/d      m*n*3
c          low clouds (cc)
c   aerosol optical thickness (taual)                n/d      m*ndim*np
c   cosecant of the solar zenith angle (csm)         n/d      m*n
c   near ir surface albedo for beam                fraction   m*ndim
c                radiation (rsirbm)
c   near ir surface albedo for diffuse             fraction   m*ndim
c                radiation (rsirdf)
c
c----- output (updated) parameters:
c
c   all-sky flux (downward-upward) (flx)           fraction   m*ndim*(np+1)
c   clear-sky flux (downward-upward) (flc)         fraction   m*ndim*(np+1)
c   all-sky direct downward ir flux at
c          the surface (fdirir)                    fraction   m*ndim
c   all-sky diffuse downward ir flux at
c          the surface (fdifir)                    fraction   m*ndim
c
c----- note: the following parameters must be specified by users:
c   aerosol single scattering albedo (ssaal)         n/d      nband
c   aerosol asymmetry factor (asyal)                 n/d      nband
c
c*************************************************************************

      implicit none

c-----Explicit Inline Directives 
 
#ifdef CRAY
#ifdef f77 
cfpp$ expand (deledd)
cfpp$ expand (sagpol)
cfpp$ expand (expmn) 
#endif 
#endif 
      _RL expmn

c-----input parameters

      integer m,n,ndim,np,ict,icb
      _RL  taucld(m,ndim,np,2),reff(m,ndim,np,2),fcld(m,ndim,np)
      _RL  tauclb(m,n,np),tauclf(m,n,np),cc(m,n,3)
      _RL  rsirbm(m,ndim),rsirdf(m,ndim)
      _RL  wh(m,n,np),taual(m,ndim,np),csm(m,n)

c-----output (updated) parameters

      _RL  flx(m,ndim,np+1),flc(m,ndim,np+1)
      _RL  fdirir(m,ndim),fdifir(m,ndim)

c-----static parameters

      integer nk,nband
      parameter (nk=10,nband=3)
      _RL  xk(nk),hk(nband,nk),ssaal(nband),asyal(nband)
      _RL  aia(nband,3),awa(nband,3),aig(nband,3),awg(nband,3)

c-----temporary array

      integer ib,ik,i,j,k
      _RL  ssacl(m,n,np),asycl(m,n,np)
      _RL  rr(m,n,np+1,2),tt(m,n,np+1,2),td(m,n,np+1,2),
     *       rs(m,n,np+1,2),ts(m,n,np+1,2)
      _RL  fall(m,n,np+1),fclr(m,n,np+1)
      _RL  fsdir(m,n),fsdif(m,n)

      _RL  tauwv,tausto,ssatau,asysto,tauto,ssato,asyto
      _RL  taux,reff1,reff2,w1,w2,g1,g2
      _RL  ssaclt(m,n),asyclt(m,n)
      _RL  rr1t(m,n),tt1t(m,n),td1t(m,n),rs1t(m,n),ts1t(m,n)
      _RL  rr2t(m,n),tt2t(m,n),td2t(m,n),rs2t(m,n),ts2t(m,n)

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

      data xk/            
     1  0.0010, 0.0133, 0.0422, 0.1334, 0.4217,            
     2  1.334,  5.623,  31.62,  177.8,  1000.0/  

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

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

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

      data ssaal,asyal/nband*0.999,nband*0.850/
 
c-----coefficients for computing the single scattering albedo of
c     ice clouds from ssa=1-(aia(*,1)+aia(*,2)*reff+aia(*,3)*reff**2)

      data aia/
     1  .08938331, .00215346,-.00000260,
     2  .00299387, .00073709, .00000746,
     3 -.00001038,-.00000134, .00000000/

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

      data awa/
     1  .01209318,-.00019934, .00000007,
     2  .01784739, .00088757, .00000845,
     3 -.00036910,-.00000650,-.00000004/

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

      data aig/
     1  .84090400, .76098937, .74935228,
     2  .00126222, .00141864, .00119715,
     3 -.00000385,-.00000396,-.00000367/

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

      data awg/
     1  .83530748, .74513197, .79375035,
     2  .00257181, .01370071, .00832441,
     3  .00005519,-.00038203,-.00023263/

c-----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)=rsirbm(i,j)
         rs(i,j,np+1,1)=rsirdf(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
       enddo
      enddo

c-----integration over spectral bands

      do 100 ib=1,nband

c-----compute cloud single scattering albedo and asymmetry factor
c     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=taucld(i,j,k,1)+taucld(i,j,k,2)
          if (taux.gt.0.05 .and. fcld(i,j,k).gt.0.01) then

           reff1=min(reff(i,j,k,1),130. _d 0)
           reff2=min(reff(i,j,k,2),20.0 _d 0)

           w1=(1.-(aia(ib,1)+(aia(ib,2)+
     *         aia(ib,3)*reff1)*reff1))*taucld(i,j,k,1)
           w2=(1.-(awa(ib,1)+(awa(ib,2)+
     *         awa(ib,3)*reff2)*reff2))*taucld(i,j,k,2)
           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)

          endif

         enddo
        enddo

        do j=1,n
         do i=1,m
            ssacl(i,j,k)=ssaclt(i,j)
         enddo
        enddo
        do j=1,n
         do i=1,m
            asycl(i,j,k)=asyclt(i,j)
         enddo
        enddo

       enddo

c-----integration over the k-distribution function

         do 200 ik=1,nk

          do 300 k= 1, np

           do j= 1, n
            do i= 1, m

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

c-----compute total optical thickness, single scattering albedo,
c     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)
 
c-----compute reflectance and transmittance for cloudless layers

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

c            if (ssato .gt. 0.001) then

c             ssato=min(ssato,0.999999 _d 0)
c             asyto=asysto/(ssato*tauto)

c             call deledd(tauto,ssato,asyto,csm(i,j), 
c     *                   rr1t(i,j),tt1t(i,j),td1t(i,j))

c             call sagpol (tauto,ssato,asyto,rs1t(i,j),ts1t(i,j))

c            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

c            endif

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

             if (tauclb(i,j,k) .lt. 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 _d 0)
              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 _d 0)
              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))

             endif

            enddo
           enddo

c-----the following code appears to be awkward, but it is efficient
c     in certain parallel processors.

           do j=1,n
            do i=1,m
               rr(i,j,k,1)=rr1t(i,j)
            enddo
           enddo
           do j=1,n
            do i=1,m
               tt(i,j,k,1)=tt1t(i,j)
            enddo
           enddo
           do j=1,n
            do i=1,m
               td(i,j,k,1)=td1t(i,j)
            enddo
           enddo
           do j=1,n
            do i=1,m
               rs(i,j,k,1)=rs1t(i,j)
            enddo
           enddo
           do j=1,n
            do i=1,m
               ts(i,j,k,1)=ts1t(i,j)
            enddo
           enddo
 
           do j=1,n
            do i=1,m
               rr(i,j,k,2)=rr2t(i,j)
            enddo
           enddo
           do j=1,n
            do i=1,m
               tt(i,j,k,2)=tt2t(i,j)
            enddo
           enddo
           do j=1,n
            do i=1,m
               td(i,j,k,2)=td2t(i,j)
            enddo
           enddo
           do j=1,n
            do i=1,m
               rs(i,j,k,2)=rs2t(i,j)
            enddo
           enddo
           do j=1,n
            do i=1,m
               ts(i,j,k,2)=ts2t(i,j)
            enddo
           enddo

 300  continue

c-----flux calculations
 
       do k= 1, np+1
        do j= 1, n
         do i= 1, m
          fclr(i,j,k) = 0.
          fall(i,j,k) = 0.
         enddo
        enddo
       enddo
       do j= 1, n
        do i= 1, m
         fsdir(i,j) = 0.
         fsdif(i,j) = 0.
        enddo
       enddo

        call cldflx (m,n,np,ict,icb,cc,rr,tt,td,rs,ts,
     *               fclr,fall,fsdir,fsdif)

       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)
         enddo
        enddo
        do j= 1, n
         do i= 1, m
          flc(i,j,k) = flc(i,j,k)+fclr(i,j,k)*hk(ib,ik)
         enddo
        enddo
       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)
        enddo
       enddo

  200 continue
  100 continue
 
      return
      end

c************************************************************************

      subroutine soluv (m,n,ndim,np,oh,dp,taucld,tauclb,tauclf,reff,
     *  ict,icb,fcld,cc,taual,csm,rsuvbm,rsuvdf,flx,flc
     * ,fdirpar,fdifpar,fdiruv,fdifuv)

c************************************************************************
c  compute solar fluxes in the uv+visible region. the spectrum is
c  grouped into 8 bands:
c  
c              Band     Micrometer
c
c       UV-C    1.     .175 - .225
c               2.     .225 - .245
c                      .260 - .280
c               3.     .245 - .260
c
c       UV-B    4.     .280 - .295
c               5.     .295 - .310
c               6.     .310 - .320
c      
c       UV-A    7.     .320 - .400
c      
c       PAR     8.     .400 - .700
c
c----- Input parameters:                            units      size
c
c   number of soundings in zonal direction (m)       n/d        1
c   number of soundings in meridional direction (n)  n/d        1
c   maximum number of soundings in                   n/d        1
c           meridional direction (ndim)
c   number of atmospheric layers (np)                n/d        1
c   layer ozone content (oh)                      (cm-atm)stp m*n*np
c   layer pressure thickness (dp)                    mb       m*n*np
c   cloud optical thickness (taucld)                 n/d      m*ndim*np*2
c          index 1 for ice paticles
c          index 2 for liquid particles
c   scaled cloud optical thickness                   n/d      m*n*np
c          for beam radiation (tauclb)
c   scaled cloud optical thickness                   n/d      m*n*np
c          for diffuse radiation  (tauclf)
c   effective cloud-particle size (reff)           micrometer m*ndim*np*2
c          index 1 for ice paticles
c          index 2 for liquid particles
c   level indiex separating high and                 n/d      m*n
c          middle clouds (ict)
c   level indiex separating middle and               n/d      m*n
c          low clouds (icb)
c   cloud amount (fcld)                            fraction   m*ndim*np
c   cloud amounts of high, middle, and               n/d      m*n*3
c          low clouds (cc)
c   aerosol optical thickness (taual)                n/d      m*ndim*np
c   cosecant of the solar zenith angle (csm)         n/d      m*n
c   uv+par surface albedo for beam                 fraction   m*ndim
c           radiation (rsuvbm)
c   uv+par surface albedo for diffuse              fraction   m*ndim
c           radiation (rsuvdf)
c
c----- output (updated) parameters:
c
c   all-sky net downward flux (flx)                fraction   m*ndim*(np+1)
c   clear-sky net downward flux (flc)              fraction   m*ndim*(np+1)
c   all-sky direct downward par flux at
c          the surface (fdirpar)                   fraction   m*ndim
c   all-sky diffuse downward par flux at
c          the surface (fdifpar)                   fraction   m*ndim
c
c----- note: the following parameters must be specified by users:
c
c   aerosol single scattering albedo (ssaal)         n/d        1
c   aerosol asymmetry factor (asyal)                 n/d        1
c
*
c***********************************************************************

      implicit none

c-----Explicit Inline Directives  
  
#ifdef CRAY 
#ifdef f77  
cfpp$ expand (deledd) 
cfpp$ expand (sagpol) 
#endif  
#endif

c-----input parameters

      integer m,n,ndim,np,ict,icb
      _RL  taucld(m,ndim,np,2),reff(m,ndim,np,2),fcld(m,ndim,np)
      _RL  tauclb(m,n,np),tauclf(m,n,np),cc(m,n,3)
      _RL  oh(m,n,np),dp(m,n,np),taual(m,ndim,np)
      _RL  rsuvbm(m,ndim),rsuvdf(m,ndim),csm(m,n)

c-----output (updated) parameter

      _RL  flx(m,ndim,np+1),flc(m,ndim,np+1)
      _RL  fdirpar(m,ndim),fdifpar(m,ndim)
      _RL  fdiruv(m,ndim),fdifuv(m,ndim)

c-----static parameters

      integer nband
      parameter (nband=8)
      _RL  hk(nband),xk(nband),ry(nband)
      _RL  asyal(nband),ssaal(nband),aig(3),awg(3)

c-----temporary array

      integer i,j,k,ib
      _RL  taurs,tauoz,tausto,ssatau,asysto,tauto,ssato,asyto
      _RL  taux,reff1,reff2,g1,g2,asycl(m,n,np)
      _RL  td(m,n,np+1,2),rr(m,n,np+1,2),tt(m,n,np+1,2),
     *       rs(m,n,np+1,2),ts(m,n,np+1,2)
      _RL  fall(m,n,np+1),fclr(m,n,np+1),fsdir(m,n),fsdif(m,n)
      _RL  asyclt(m,n)
      _RL  rr1t(m,n),tt1t(m,n),td1t(m,n),rs1t(m,n),ts1t(m,n)
      _RL  rr2t(m,n),tt2t(m,n),td2t(m,n),rs2t(m,n),ts2t(m,n)

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

      data hk/.00057, .00367, .00083, .00417,
     *        .00600, .00556, .05913, .39081/

c-----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/

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

      data ry /.00604, .00170, .00222, .00132,
     *         .00107, .00091, .00055, .00012/

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

      data ssaal,asyal/nband*0.999,nband*0.850/

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

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

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

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

c-----initialize surface reflectances and transmittances
            
      do j= 1, n
       do i= 1, m                    
         rr(i,j,np+1,1)=rsuvbm(i,j)
         rr(i,j,np+1,2)=rsuvbm(i,j)
         rs(i,j,np+1,1)=rsuvdf(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
       enddo
      enddo

c-----compute cloud asymmetry factor for a mixture of
c     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=taucld(i,j,k,1)+taucld(i,j,k,2)
          if (taux.gt.0.05 .and. fcld(i,j,k).gt.0.01) then

           reff1=min(reff(i,j,k,1),130. _d 0)
           reff2=min(reff(i,j,k,2),20.0 _d 0)

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

          endif

        enddo
       enddo

       do j=1,n
        do i=1,m
           asycl(i,j,k)=asyclt(i,j)
        enddo
       enddo

      enddo
 
      do j=1,n
       do i=1,m
        fdiruv(i,j)=0.0
        fdifuv(i,j)=0.0
       enddo
      enddo
           
c-----integration over spectral bands

      do 100 ib=1,nband

       do 300 k= 1, np

        do j= 1, n
         do i= 1, m

c-----compute ozone and rayleigh optical thicknesses

          taurs=ry(ib)*dp(i,j,k)
          tauoz=xk(ib)*oh(i,j,k)
 
c-----compute total optical thickness, single scattering albedo,
c     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)

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

          tauto=tausto
          ssato=ssatau/tauto+1.0e-8
          ssato=min(ssato,0.999999 _d 0)
          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))

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

          if (tauclb(i,j,k) .lt. 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 _d 0)
           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 _d 0)
           asyto=(asysto+asycl(i,j,k)*tauclf(i,j,k))/(ssato*tauto)

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

          endif

         enddo
        enddo

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

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

 300  continue

c-----flux calculations
 
       do k= 1, np+1
        do j= 1, n
         do i= 1, m
          fclr(i,j,k) = 0.
          fall(i,j,k) = 0.
         enddo
        enddo
       enddo
       do j= 1, n
        do i= 1, m
         fsdir(i,j) = 0.
         fsdif(i,j) = 0.
        enddo
       enddo
        call cldflx (m,n,np,ict,icb,cc,rr,tt,td,rs,ts,
     *               fclr,fall,fsdir,fsdif)

       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)
         enddo
        enddo
        do j= 1, n
         do i= 1, m
          flc(i,j,k)=flc(i,j,k)+fclr(i,j,k)*hk(ib)
         enddo
        enddo
       enddo

       if(ib.eq.nband) 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)
         enddo
        enddo
       else
         do j=1,n
          do i=1,m
           fdiruv(i,j)=fdiruv(i,j)+fsdir(i,j)*hk(ib)
           fdifuv(i,j)=fdifuv(i,j)+fsdif(i,j)*hk(ib)
         enddo
        enddo
       endif

 100  continue

      return
      end

c*********************************************************************

      subroutine deledd(tau,ssc,g0,csm,rr,tt,td)

c*********************************************************************
c
c-----uses the delta-eddington approximation to compute the
c     bulk scattering properties of a single layer
c     coded following King and Harshvardhan (JAS, 1986)
c
c  inputs:
c
c     tau: the effective optical thickness
c     ssc: the effective single scattering albedo
c     g0:  the effective asymmetry factor
c     csm: the effective secant of the zenith angle
c
c  outputs:
c
c     rr: the layer reflection of the direct beam
c     tt: the layer diffuse transmission of the direct beam
c     td: the layer direct transmission of the direct beam

c*********************************************************************

      implicit none

c-----Explicit Inline Directives  
  
#ifdef CRAY 
#ifdef f77  
cfpp$ expand (expmn)
#endif  
#endif
      _RL expmn

      _RL  zero,one,two,three,four,fourth,seven,tumin
      parameter (one=1., three=3.)
      parameter (seven=7., two=2.)
      parameter (four=4., fourth=.25)
      parameter (zero=0., tumin=1.e-20)

c-----input parameters
      _RL  tau,ssc,g0,csm

c-----output parameters
      _RL  rr,tt,td

c-----temporary parameters

      _RL  zth,ff,xx,taup,sscp,gp,gm1,gm2,gm3,akk,alf1,alf2,
     *     all1,bll,st7,st8,cll,dll,fll,ell,st1,st2,st3,st4
c
                zth = one / csm
 
c  delta-eddington scaling of single scattering albedo,
c  optical thickness, and asymmetry factor,
c  K & H eqs(27-29)

                ff  = g0*g0
                xx  = one-ff*ssc
                taup= tau*xx
                sscp= ssc*(one-ff)/xx
                gp  = g0/(one+g0)
 
c  extinction of the direct beam transmission
 
                td  = expmn(-taup*csm)
 
c  gamma1, gamma2, gamma3 and gamma4. see table 2 and eq(26) K & H
c  ssc and gp are the d-s single scattering
c  albedo and asymmetry factor.

                xx  =  three*gp 
                gm1 =  (seven - sscp*(four+xx))*fourth
                gm2 = -(one   - sscp*(four-xx))*fourth
                gm3 =  (two   -        zth*xx )*fourth
 
c  akk is k as defined in eq(25) of K & H
 
                xx  = gm1-gm2
                akk = sqrt((gm1+gm2)*xx)
 
c   alf1 and alf2 are alpha1 and alpha2 from eqs (23) & (24) of K & H
 
                alf1 = gm1 - gm3 * xx
                alf2 = gm2 + gm3 * xx
 
c  all1 is last term in eq(21) of K & H
c  bll is last term in eq(22) of K & H
 
                xx  = akk * two
                all1 = (gm3 - alf2 * zth    )*xx*td
                bll = (one - gm3 + alf1*zth)*xx
 
                xx  = akk * zth
                st7 = one - xx
                st8 = one + xx
 
                xx  = akk * gm3
                cll = (alf2 + xx) * st7
                dll = (alf2 - xx) * st8
 
                xx  = akk * (one-gm3)
                fll = (alf1 + xx) * st8
                ell = (alf1 - xx) * st7
 
                st3 = max(abs(st7*st8),tumin)
                st3 = sign (st3,st7)
 
                st2 = expmn(-akk*taup)
                st4 = st2 * st2
 
                st1 =  sscp / ((akk+gm1 + (akk-gm1)*st4) * st3)
 
c   rr is r-hat of eq(21) of K & H
c   tt is diffuse part of t-hat of eq(22) of K & H
 
                rr =   ( cll-dll*st4    -all1*st2)*st1
                tt = - ((fll-ell*st4)*td-bll*st2)*st1
 
                rr = max(rr,zero)
                tt = max(tt,zero)

      return
      end

c*********************************************************************

      subroutine sagpol(tau,ssc,g0,rll,tll)

c*********************************************************************
c-----transmittance (tll) and reflectance (rll) of diffuse radiation
c     follows Sagan and Pollock (JGR, 1967).
c     also, eq.(31) of Lacis and Hansen (JAS, 1974).
c
c-----input parameters:
c
c      tau: the effective optical thickness
c      ssc: the effective single scattering albedo
c      g0:  the effective asymmetry factor
c
c-----output parameters:
c
c      rll: the layer reflection of diffuse radiation
c      tll: the layer transmission of diffuse radiation
c
c*********************************************************************

      implicit none

c-----Explicit Inline Directives  
  
#ifdef CRAY 
#ifdef f77  
cfpp$ expand (expmn)
#endif  
#endif
      _RL expmn

      _RL  one,three,four
      parameter (one=1., three=3., four=4.)

c-----output parameters:

      _RL  tau,ssc,g0

c-----output parameters:

      _RL  rll,tll

c-----temporary arrays

      _RL  xx,uuu,ttt,emt,up1,um1,st1

             xx  = one-ssc*g0
             uuu = sqrt( xx/(one-ssc))
             ttt = sqrt( xx*(one-ssc)*three )*tau
             emt = expmn(-ttt)
             up1 = uuu + one
             um1 = uuu - one
             xx  = um1*emt
             st1 = one / ((up1+xx) * (up1-xx))
             rll = up1*um1*(one-emt*emt)*st1
             tll = uuu*four*emt         *st1

      return
      end

c*******************************************************************

      function expmn(fin)

c*******************************************************************
c compute exponential for arguments in the range 0> fin > -10.
c*******************************************************************
      implicit none
      _RL  fin,expmn

      _RL one,expmin,e1,e2,e3,e4
      parameter (one=1.0, expmin=-10.0)
      parameter (e1=1.0,        e2=-2.507213e-1)
      parameter (e3=2.92732e-2, e4=-3.827800e-3)

      if (fin .lt. expmin) fin = expmin
      expmn = ((e4*fin + e3)*fin+e2)*fin+e1
      expmn = expmn * expmn
      expmn = one / (expmn * expmn)

      return
      end

c*******************************************************************

      subroutine cldflx (m,n,np,ict,icb,cc,rr,tt,td,rs,ts,
     *           fclr,fall,fsdir,fsdif)

c*******************************************************************
c  compute upward and downward fluxes using a two-stream adding method
c  following equations (3)-(5) of Chou (1992, JAS).
c
c  clouds are grouped into high, middle, and low clouds which are 
c  assumed randomly overlapped. It involves eight sets of calculations.
c  In each set of calculations, each atmospheric layer is homogeneous,
c  either with clouds or without clouds.

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

      implicit none

c-----input parameters

      integer m,n,np,ict,icb

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

c-----temporary array

      integer i,j,k,ih,im,is
      _RL  rra(m,n,np+1,2,2),tta(m,n,np+1,2,2),tda(m,n,np+1,2,2)
      _RL  rsa(m,n,np+1,2,2),rxa(m,n,np+1,2,2)
      _RL  ch(m,n),cm(m,n),ct(m,n),flxdn(m,n,np+1)
      _RL  fdndir(m,n),fdndif(m,n),fupdif
      _RL  denm,xx

c-----output parameters

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

c-----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
        enddo
       enddo
      enddo

       do j=1,n
        do i=1,m
           fsdir(i,j)=0.0
           fsdif(i,j)=0.0
        enddo
       enddo

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

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

      do 10 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)
        enddo
       enddo

       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)
         enddo
        enddo
       enddo

c-----for middle clouds

      do 10 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
         enddo
        enddo
       enddo

  10  continue

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

      do 20 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)
        enddo
       enddo

       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)
         enddo
        enddo
       enddo

c-----for middle clouds

      do 20 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
         enddo
        enddo
       enddo

  20  continue

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

      do 100 ih=1,2

c-----clear portion 

         if(ih.eq.1) then
           do j=1,n
            do i=1,m
             ch(i,j)=1.0-cc(i,j,1)
            enddo
           enddo

          else

c-----cloudy portion

           do j=1,n
            do i=1,m
             ch(i,j)=cc(i,j,1)
            enddo
           enddo

          endif

      do 100 im=1,2

c-----clear portion

         if(im.eq.1) then

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

         else

c-----cloudy portion

           do j=1,n
            do i=1,m
              cm(i,j)=ch(i,j)*cc(i,j,2) 
            enddo
           enddo

         endif

      do 100 is=1,2

c-----clear portion

         if(is.eq.1) then

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

         else

c-----cloudy portion

           do j=1,n
            do i=1,m
             ct(i,j)=cm(i,j)*cc(i,j,3)
            enddo
           enddo

         endif

c-----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
         enddo
        enddo
       enddo

c-----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
         enddo
        enddo
       enddo

c-----compute fluxes following eq (5) of Chou (1992)
 
c     fdndir is the direct  downward flux
c     fdndif is the diffuse downward flux
c     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
        enddo
       enddo
      enddo

       do j=1, n
        do i=1, m
         flxdn(i,j,1)=1.0-rra(i,j,1,im,is)
        enddo
       enddo

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

       do k=1,np+1
        do j=1,n
         do i=1,m
           if(ih.eq.1 .and. im.eq.1 .and. is.eq.1) then
             fclr(i,j,k)=flxdn(i,j,k)
           endif
             fall(i,j,k)=fall(i,j,k)+flxdn(i,j,k)*ct(i,j)
         enddo
        enddo
       enddo

        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)
         enddo
        enddo

  100 continue

      return
      end

c*****************************************************************

      subroutine flxco2(m,n,np,swc,swh,csm,df)

c*****************************************************************

c-----compute the reduction of clear-sky downward solar flux
c     due to co2 absorption.

      implicit none

c-----input parameters

      integer m,n,np
      _RL  csm(m,n),swc(m,n,np+1),swh(m,n,np+1),cah(22,19)

c-----output (undated) parameter

      _RL  df(m,n,np+1)

c-----temporary array

      integer i,j,k,ic,iw 
      _RL  xx,clog1,wlog,dc,dw,x1,x2,y2

c********************************************************************
c-----include co2 look-up table

#include "cah-dat.h"
c     save cah
 
c********************************************************************
c-----table look-up for the reduction of clear-sky solar
c     radiation due to co2. The fraction 0.0343 is the
c     extraterrestrial solar flux in the co2 bands.

      do k= 2, np+1
       do j= 1, n
        do i= 1, m
          xx=1./.3
          clog1=log10(swc(i,j,k)*csm(i,j))
          wlog=log10(swh(i,j,k)*csm(i,j))
          ic=int( (clog1+3.15)*xx+1.)
          iw=int( (wlog+4.15)*xx+1.)
          if(ic.lt.2)ic=2
          if(iw.lt.2)iw=2
          if(ic.gt.22)ic=22
          if(iw.gt.19)iw=19     
          dc=clog1-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.lt.y2) x1=y2
          df(i,j,k)=df(i,j,k)+0.0343*(x1-y2)
        enddo     
       enddo      
      enddo      

      return
      end