subroutine radsw
c************
c a simplified solar radiation scheme derived from Fu-Liou scheme (1993).
c             FSW=S0*Ca(1)(1+bs(1))
c             FSWds=S0*Ca(2)(1+bs(2))
c      where
c             Ca(i=1,2)=a(i=1,2)*cosZ+b(i=1,2)
c             bs(i=1,2)=c(i=1,2)*ALBDs
c
c     cosZ    : cosine solar zenith angle
c     FSWds   : downward solar radiation at the surface
c     FSWus   : solar reflection at the surface
c     FSWut   : upward solar radiation at the top of the atmosphere
c     SolarC  : solar constant
c     S0      : S0 = SolarC * cosZ
c     ALBDs   : surface albedo
c     Ca(1)   : atmoshperic column absorption including clouds
c     Ca(2)   : ratio of surface solar irradiance and S0
c     bs(1)   : first order contribution due to surface reflection
c     bs(2)   : first order contribution due to surface reflection
c
c************
c Solar radiation at TOA (S0), without diurnal cycle
c Astronomical formulae used:
c 1.  solar declination angle as a function of dayofyear
c 2.  correction(<3.4%) to solar constant due to the
c     ellipticity of the earth orbit, function of dayofyear  Dickinson(1983)
c
c --Local Variables & Functions:
c    SolarC : solar constant
c    cosZ   :  zenith angle in cosine, cosZ=0 at night
c    cosZ1  : day average of cosZ
c    cosZ2  : day average of cosZ**2
c    lamda  : latitude
c    delta  : solar inclination
c    H      : hour angle at sunset and sunrise
c
c --Altered Common Variables:
c FSW          !total SW absorbed by column   W/m**2
c
      implicit none
c
c Horizontal resolution parameters
c
      integer NX,NY,NXY
      parameter(NX=64,NY=33)   !domain is periodic in x, to 60 N/S
      parameter(NXY=NX*NY)
c
      integer clt
      parameter(clt=3)  !cloud types
      real*8 cosZ                             !solar zenith angle
     1      ,cld(NX,NY,0:clt)           !cloud cover
c
      real*8 S0(NX,NY)                      !incoming solar at top
     1       ,FSWds(NX,NY),FSWus(NX,NY)      !down/up SW fluxes at surface
     2       ,FSWut(NX,NY),FSW(NX,NY)        !SW at top; SW absorbed by atmo.
             ,dayofyear
c
      integer i,j,k
      real*8 cosZ1,cosZ2,SolarC
      real*8 lamda,sinlamda,sindelta,coslamda
     2,      cosdelta,sinH,sin2H,cosH
     3,      dtheta,delta,H,ds
      real*8 dayspring,dayperigee,ecce,theta0max
      real*8 Ca(0:clt,2),bs(0:clt,2),a(0:clt,2),b(0:clt,2),c(0:clt,2)

      data (a(i,1),i=0,clt)/
     *              -0.6991568D-01,  0.2862155D-01, -0.3275377D-01,
     *              -0.3075491D-01/
      data (b(i,1),i=0,clt)/
     *               0.2489715D+00,  0.1758697D+00,  0.2177709D+00,
     *               0.2603912D+00/
      data (a(i,2),i=0,clt)/
     *               0.1344528D+00,  0.1476796D+00,  0.2950021D+00,
     *               0.2853160D+00/
      data (b(i,2),i=0,clt)/
     *               0.6461399D+00,  0.8235896D-01,  0.3858907D+00,
     *               0.1503189D+00/
      data (c(i,1),i=0,clt)/
     *               0.1659058D+00,  0.6383526D-01,  0.1763771D+00,
     *               0.6906892D-01/
      data (c(i,2),i=0,clt)/
     *               0.7897843D-01,  0.7440951D+00,  0.2550164D+00,
     *               0.5797699D+00/
c*****************************************************************


c- orbital parameters
      theta0max=23.447                 !degree; max solar dec. angle
      dayspring = 81.                  !Julian day of Vernal Equinox
      dayperigee = 3.                  !day of year closest to the sun
      ecce = 0.034                     !related to eccentricity of earth orbit

      sindelta=dsin(theta0max/180.*pi) *                 !Zeng (1994)
     &       dsin(2.*pi*(dayofyear-dayspring)/360.)    !1yr=360day for QTCM
      delta=dasin(sindelta)
      cosdelta=dcos(delta)
      SolarC =1370.*(one+ecce*                 ! Dickinson(1983)
     &              dcos(2.*pi*(dayofyear-dayperigee)/365.))
      dtheta=dy/Rearth
      do j=1,NY
        lamda=(j-(NY+1)/2)*dtheta
        sinlamda=dsin(lamda)
        coslamda=dsqrt(one-sinlamda**2)
        cosH=-sinlamda*sindelta/(coslamda*cosdelta)
        cosH=dmax1(cosH,-one)       !polar; no night; halfday length=pi
        cosH=dmin1(cosH,one)        !polar; no day; halfday length=0
        H=dacos(cosH)
        sinH=dsqrt(one-cosH**2)
        sin2H=2*sinH*cosH
        cosZ1=(sinlamda*sindelta*H
     2             +coslamda*cosdelta*sinH)/pi
        cosZ2=((sinlamda*sindelta)**2*H
     2           +2*sinlamda*coslamda*sindelta*cosdelta*sinH
     3           +(coslamda*cosdelta)**2*(H/2.+sin2H/4.))/pi
        do k=0,clt
           Ca(k,1)=a(k,1)*cosZ2+b(k,1)*cosZ1
           Ca(k,2)=a(k,2)*cosZ2+b(k,2)*cosZ1
        enddo
      do i=1,NX
c******************************* empirical formula from Fu-Liou scheme
        do k=0,clt
           bs(k,1)=c(k,1)*ALBDs(i,j)
           bs(k,2)=c(k,2)*ALBDs(i,j)
        enddo
cc
        FSW(i,j)=0
        FSWds(i,j)=0
        do k=0,clt
           FSW(i,j)=FSW(i,j)+cld(i,j,k)*Ca(k,1)*(1.+bs(k,1))
           FSWds(i,j)=FSWds(i,j)+cld(i,j,k)*Ca(k,2)*(1.+bs(k,2))
        enddo
        S0(i,j)=SolarC*cosZ1
        FSW(i,j)=SolarC*FSW(i,j)                    !atmospheric absorption
        FSW(i,j)=FSW(i,j)                           ! sensitivity test
        FSWds(i,j)=SolarC*FSWds(i,j)                ! surface solar irradiance
        FSWus(i,j)=FSWds(i,j)*ALBDs(i,j)
        FSWut(i,j)=S0(i,j)-FSWds(i,j)-FSW(i,j)+FSWus(i,j)
      end do
      end do

      return
      end