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     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    cosZ   :  zenith angle in cosine, cosZ=0 at night
c    lamda  : latitude
c    delta  : solar inclination
c    hh     : hour angle, h=0 at noon
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
c --Altered Common Variables:
c FSW          !total SW absorbed by column   W/m**2

c --Local Variables & Functions:
      integer i,j,k,m,maxb
      parameter(maxb=72)
      real*8 lamda,delta,hh
      real*8 sinlamda(ny),sindelta,coslamda(ny)
     2,      cosdelta,coshh(nx,maxb),SolarC
      real*8 dayspring,dayperigee,ecce,theta0max,dtheta
      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

      if(timeofday .eq. 0.)then
        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
        sindelta=dsin(delta)
        cosdelta=dsqrt(one-sindelta**2)
        if(dayofyear .eq. 1)then
          if((86400/dt) .gt. maxb)
     2            write(*,*)'hour angle array is too small, change maxb'
          do j=1,ny
            lamda=(j-(NY+1)/2)*dtheta
            sinlamda(j)=dsin(lamda)
            coslamda(j)=dsqrt(one-sinlamda(j)**2)
          enddo
          m=86400./dt
          do k=1,m
          do i=1,NX
            hh=2*pi*k/m-pi+(i-1/2.)*dx/Rearth
            coshh(i,k)=dcos(hh)
          enddo
          enddo
        endif
      endif
      m=timeofday/dt+1
      do j=1,NY
      do i=1,NX
        cosZ=sinlamda(j)*sindelta+coslamda(j)*cosdelta*coshh(i,m)
        S0(i,j)=0.
        FSWds(i,j)=0.                              ! at night
        FSW(i,j)=0.
        FSWus(i,j)=0.
        FSWut(i,j)=0.
        if(cosZ .gt. 1.D-8)then               ! at day
c******************************* empirical formula from Fu-Liou scheme
        do k=0,clt
           Ca(k,1)=a(k,1)*cosZ+b(k,1)
           Ca(k,2)=a(k,2)*cosZ+b(k,2)
        enddo
        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*cosZ
        FSW(i,j)=S0(i,j)*FSW(i,j)                    !atmospheric absorption
        FSWds(i,j)=S0(i,j)*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)
        endif
      end do
      end do

      return
      end