SUBROUTINE radfs(ix,ix2,qs,pp,ppi,qqh2o,tt,o3qo3,tsfc,slmsk,albedo,xlat,  &
     camt,ityp,ktop,kbtm,nclds,emcld,rrcl,ttcl,coszro,taudar,ibeg,ipts,  &
     ko3,kalb,slmrf,slymrf,itimsw,itimlw,jd,r1,gmt,swh,hlw,  &
     flwup,fswup,fswdn,fswdns,fswups,flwdns,flwups)
  !
  ! $Author: CPTEC $
  ! $Date: 2000/02/24 13:15:43 $
  ! $Revision: 1.0 $
  !
  !the internal drive for gfdl radiation                       
  !this subroutine was from y.h and k.a.c (1993)               
  !and modified by q. zhao for use in the eta model            
  !                 nov. 18,  1993                              
  !                                                  
  !update: this subroutine was modified to use cloud fraction    
  !        on each model layer.                                  
  !                               qingyun  zhao   95-3-22        
  !                                                              
  !update: r1 has been added to the inputs from radtn to         
  !        compute the variation of solar constant at the top    
  !        of atmosphere with julian day in a year.              
  !                               qingyun  zhao   96-7-23        
  !   ix:  horizontal real, dimension :: of one column (=2*im-1)
  !   ix2:  =ix
  !   qs(ix):  the surface pressure (pa)
  !   pp(ix,l):  the midlayer pressures (pa)  (l is the vert. dimen.)
  !   ppi(ix,lp1) the interface pressures (pa)
  !   qqh2o(ix,l): the midlayer water vapor mixing ratio (kg/kg)
  !   tt(ix,l):  the midlayer temperature (k)
  !   o3qo3(ix,l): the midlayer ozone mixing ratio
  !   tsfc(ix):  the skin temp. (k); negative over water
  !   slmsk(ix):  the sea mask (land=0,sea=1)
  !   albedo(ix): the surface albedo (expressed as a fraction)
  !   xlat(ix):  the geodetic latitudes of each column in degrees
  !    (n.h.> 0)
  ! the following are cloud information for each cloud layer
  !                      layer=1:surface
  !                      layer=2:first layer above ground, and so on
  !   camt(ix,lp1):      cloud fraction of each cloud layer
  !   ityp(ix,lp1):      cloud type(=1: stratiform, =2:convective)
  !   ktop(ix,lp1):      height of cloud top of each cloud layer (in eta level)
  !   kbtm(ix,lp1):      bottom of each cloud layer
  !   nclds(ix):         number of cloud layers
  !   emcld(ix,lp1):     cloud emissivity
  !   rrcl(ix,nb,lp1)    cloud reflecttances for sw spectral bands
  !   ttcl(ix,nb,lp1)    cloud transmitances for sw spectral bands
  ! the above are cloud information for each cloud layer
  !
  !   coszro(ix): the cosine of the solar zenith angle
  !   taudar:  =1.0
  !   ibeg:  =1
  !   ipts:  =ix
  !   ko3:  =1 ( read in the qzone data)
  !   kalb:  =0
  !   slmrf(lp1): the interface's eta (lp1=l+1)
  !   slymrf(l):  the midlayer eta
  !   itimsw:  =1/0 (shortwave calc. are desired/not desired)
  !   itimlw:  =1/0 (longwave calc. are desired/not desired)
  !the following are additional for eta model
  !   jd: julian day in a year
  !   r1: the non-real, dimension ::al sun-earth distance
  !   gmt:hour
  !generated output required by the eta model
  !    swh: atmospheric shortwave heating rates in k/s.
  !         swh is a real array real, dimension ::ed (ncol x lm).
  !    hlw: atmospheric longwave heating rates in k/s.
  !         hlw is a real array real, dimension ::ed (ncol x lm).
  !  flwup: upward longwave flux at top of the atmosphere in w/m**2.
  !         flwup is a real array real, dimension ::ed (ncol).
  !  fswup: upward shortwave flux at top of the atmosphere in w/m**2.
  !         fswup is a real array real, dimension ::ed (ncol).
  !  fswdn: downward shortwave flux at top of the atmosphere in w/m**2.
  !         fswdn is a real array real, dimension ::ed (ncol).
  ! fswdns: downward shortwave flux at the surface in w/m**2.
  !         fswdns is a real array real, dimension ::ed (ncol).
  ! fswups: upward shortwave flux at the surface in w/m**2.
  !         fswups is a real array real, dimension ::ed (ncol).
  ! flwdns: downward longwave flux at the surface in w/m**2.
  !         flwdns is a real array real, dimension ::ed (ncol).
  ! flwups: upward longwave flux at the surface in w/m**2.
  !         flwups is a real array real, dimension ::ed (ncol).
  IMPLICIT NONE

  INCLUDE "parmeta.f90"

  !generalized for plug-compatibility -                       
  !original code was cleaned-up gfdl code...k.campana mar89
  !example for mrf:                                            
  !  ko3  =0  and o3qo3=dummy array.   (gfdl climo o3 used)    
  !  kemis=0  and hi cld emis computed here (cemis=dummy input)
  !  kalb =0  and sfc albedo over open water computed below... 
  !  kcco2=0,co2 obtained from block data                      
  !         =1,co2 computed in here --- not available yet...     
  !  slmrf = interface (levels) sigma                          
  !         slymrf= layer sigma                                  
  !   updated for yutai hou sib sw radiation....kac 6 mar 92      
  !     ocean albedo for beam set to bulk sfcalb, since           
  !        cosine zenith angle effects already there(ref:payne)   
  !        slmsk = 0.                                             
  !     snow ice albedo for beam not enhanced via cosine zenith   
  !        angle either cause valu already high (we see polar    
  !        cooling if we do beam calculation)....kac 17mar92      
  !        albedo ge .5                                           
  !    updated to obtain clear sky fluxes "on the fly" for       
  !        cloud forcing diagnostics elsewhere...kac 7aug92       
  !        see ##clr lines...radfs,lwr88,fst88,spa88 .......     
  !   updated for use new cld scheme      ......yh  dec 92        
  !     input cld may be as original in 3 domain (cld,mtop,mbot)  
  !        or in a vertical array of 18 mdl layers (cldary)       
  !     iemis=0  use the org. cld emis scheme                     
  !          =1  use temp dep. cld emis scheme                    
  !   updated to compute cld layer reflecttance and transmittance 
  !     input cld emissivity and optical thickness 'emis0,tauc0'  
  !                                       ......yh feb 93         

  INCLUDE "hcon.f90"
  INCLUDE "rdparm.f90"
  INTEGER, PARAMETER :: lngth=37*l
  INCLUDE "rnddta.f90"
  INCLUDE "co2dta.f90"
  INCLUDE "tabcom.f90"

  INTEGER,INTENT(in) :: ix, ix2      
  REAL,DIMENSION(ix2),INTENT(in)  :: qs
  REAL,DIMENSION(ix,l),INTENT(in) :: pp
  REAL,DIMENSION(ix2,l),INTENT(in) :: qqh2o,tt
  REAL,DIMENSION(ix),INTENT(in) :: tsfc,slmsk,albedo,xlat,coszro,taudar
  REAL,DIMENSION(ix,lp1),INTENT(in) :: ppi,camt,emcld
  INTEGER,DIMENSION(ix),INTENT(in) :: nclds
  INTEGER,DIMENSION(ix,lp1),INTENT(in) :: ktop,kbtm 
  REAL,DIMENSION(lp1),INTENT(in) :: slmrf
  REAL,DIMENSION(l),INTENT(in) :: slymrf 
  REAL,DIMENSION(ix,nb,lp1),INTENT(inout) :: rrcl,ttcl
  REAL,DIMENSION(ix,l),INTENT(in) :: o3qo3
  REAL,DIMENSION(ix,l),INTENT(out) ::  swh,hlw
  REAL,DIMENSION(ix),INTENT(out) :: flwup,fswup,fswdn,fswdns,flwdns,flwups,fswups

  REAL, DIMENSION(ix) :: alvbr,alnbr,alvdr,alndr
  REAL, DIMENSION(ix) :: gdfvbr,gdfnbr,gdfvdr,gdfndr
  REAL, DIMENSION(imax,lp1) :: fswl,hswl,ufl,dfl
  REAL, DIMENSION(imax,lp1) :: eqcmt,press,temp,fsw,hsw,uf,df
  REAL, DIMENSION(imax,lp1,lp1) :: cldfac
  REAL, DIMENSION(imax,l) :: rh2o,qo3,heatra
  REAL, DIMENSION(imax) :: grnflx,topflx,grdflx,coszen,tauda

  !add pressure interface
  REAL, DIMENSION(nb) :: abcff,pwts
  REAL :: cfco2,cfo3,reflo3,rrayav
  COMMON /swrsav/ abcff,pwts,cfco2,cfo3,reflo3,rrayav

  !vector temporaries for cloud calc.      *
  REAL, DIMENSION(imax) :: jjrow,do3v,do3vp,tthan

  !seasonal climatologies of o3 (obtained from a previously run
  !code which interpolates o3 to user vertical coordinate).
  !defined as 5 deg lat means n.p.->s.p.
  REAL, DIMENSION(37,l) :: dduo3n,ddo3n2,ddo3n3,ddo3n4
  COMMON /savmem/dduo3n,ddo3n2,ddo3n3,ddo3n4 !winter,spring,summer,fall

  REAL, DIMENSION(lngth) :: rad1, rad2, rad3, rad4
  EQUIVALENCE (rad1(1),dduo3n(1,1)),(rad2(1),ddo3n2(1,1))
  EQUIVALENCE (rad3(1),ddo3n3(1,1)),(rad4(1),ddo3n4(1,1))

  REAL, DIMENSION(21,20) :: albd
  REAL, DIMENSION(20) :: za
  REAL, DIMENSION(21) :: trn
  REAL, DIMENSION(19) :: dza
  COMMON / ssalb / albd,za,trn,dza

  REAL, DIMENSION(imax) :: alvb,alnb,alvd,alnd,gdfvb,gdfnb,gdfvd,gdfnd,sfcalb
  REAL, DIMENSION(imax,lp1) :: xamt,mtopsw,mbtmsw
  REAL, DIMENSION(5) :: cao3sw,cah2sw,cbsw
  REAL :: emisp,emist,xlatt,xlatp,q19001,hp98,h3m6,hp75,h6m2,hp537,h74e1,h15e1
  REAL :: q14330,hp2,twenty,hnine,degrad,hsigma,daysec,rco2
  COMMON /rdfsav/ emisp,emist,xlatt,xlatp,q19001,hp98,h3m6,  &
       hp75,h6m2,hp537,h74e1,h15e1,q14330,hp2,twenty,hnine,  &
       degrad,hsigma,daysec,rco2, cao3sw,cah2sw,cbsw

  REAL :: solc,rsin1,rcos1,rcos2
  COMMON /astsav/ solc,rsin1,rcos1,rcos2

  !solc,the solar constant is scaled to a more current value.
  !i.e. if solc=2.0 ly/min then ssolar=1.96 ly/min.
  !re-computed cause ssolar overwritten as part of scratch common
  !note: xlat is in degree here
  REAL, PARAMETER :: year=365.25
  REAL, PARAMETER :: rlag=14.8125
  REAL, PARAMETER :: tpi=6.283185308
  REAL, PARAMETER :: pi=3.1415927
  REAL, PARAMETER :: sc=2.

  INTEGER :: i,ipts,ir,ibeg,iq,itimlw,itimsw,ityp,jd,jx,k,kalb,ko3,n
  REAL :: r1,ssolar,date,gmt,rang,th2,zen,dzen,alb1,alb2,albd0,alvd1,alnd1
  REAL :: rrvco2,rrco2 


  ! special note: the solar constant is reduced extra 3 percent to account
  !               for the lack of aerosols in the shortwave radiation
  !               parameterization.       q. zhao    96-7-23
  solc=sc/(r1*r1)
  ssolar=solc*hp98
  ssolar=ssolar*0.97
  date=jd+gmt/24.0
  rang=tpi*(date-rlag)/year
  rsin1=SIN(rang)
  rcos1=COS(rang)
  rcos2=COS(2.0*rang)

  DO  i=1,ipts
     ir = i + ibeg - 1
     th2=hp2*xlat(ir)
     jjrow(i)=q19001-th2
     tthan(i)=(19-jjrow(i))-th2
     !note that the nmc variables are in mks (thus pressure is in
     !centibars)while all gfdl variables are in cgs units
     sfcalb(i) = albedo(ir)
     !now put sfc temp,pressures, zenith angle into sw common block...
     !note: all pressures input from the eta model are in pa
     !      the unit for press is micro bar
     !      surface temperature are negative over oceans in the eta model
     press(i,lp1)=qs(ir)*10.0
     temp(i,lp1)=ABS(tsfc(ir))
     coszen(i) = coszro(ir)
     tauda(i) = taudar(ir)
  END DO
  !all gfdl variables have k=1 at the top of the atmosphere.nmc
  !eta model has the same structure
  DO  k=1,l
     DO  i=1,ipts
        ir = i + ibeg - 1
        !now put temp,pressures, into sw common block..........
        temp(i,k) = tt(ir,k)
        press(i,k) = 10.0 * pp(ir,k)
        !store lyr moisture and add to sw common block
        rh2o(i,k)=qqh2o(ir,k)
        IF(rh2o(i,k) < h3m6) rh2o(i,k)=h3m6
     END DO
  END DO
  IF (ko3 /= 0) THEN
     DO  k=1,l
        DO  i=1,ipts
           qo3(i,k) = o3qo3(i+ibeg-1,k)
        END DO
     END DO
  END IF
  IF (kalb <= 0) THEN
     !the following code gets albedo from payne,1972 tables if
     !1) open sea point (slmsk=1);2) kalb=0
     iq=INT(twenty*hp537+one)
     DO  i=1,ipts
        IF(coszen(i) > 0.0 .AND. slmsk(i+ibeg-1) > 0.5) THEN
           zen=degrad*ACOS(MAX(coszen(i),0.0))
           IF(zen >= h74e1) jx=INT(haf*(hninety-zen)+one)
           IF(zen < h74e1.AND.zen >= fifty) jx=INT(quartr*(h74e1-zen)+hnine)
           IF(zen < fifty) jx=INT(hp1*(fifty-zen)+h15e1)
           dzen=-(zen-za(jx))/dza(jx)
           alb1=albd(iq,jx)+dzen*(albd(iq,jx+1)-albd(iq,jx))
           alb2=albd(iq+1,jx)+dzen*(albd(iq+1,jx+1)-albd(iq+1,jx))
           sfcalb(i)=alb1+twenty*(alb2-alb1)*(hp537-trn(iq))
        END IF
     END DO
  END IF
  IF (ko3 <= 0) THEN
     !compute climatological zonal mean ozone,
     !seasonal and spatial interpolation done below.
     DO  k=1,l
        DO  i=1,ipts
           do3v(i)  = dduo3n(jjrow(i),k) + rsin1*ddo3n2(jjrow(i),k)  &
                +rcos1*ddo3n3(jjrow(i),k) +rcos2*ddo3n4(jjrow(i),k)
           do3vp(i) = dduo3n(jjrow(i)+1,k) + rsin1*ddo3n2(jjrow(i)+1,k)  &
                +rcos1*ddo3n3(jjrow(i)+1,k) +rcos2*ddo3n4(jjrow(i)+1,k)
           !now latitudinal interpolation, and
           !convert o3 into mass mixing ratio(original data mpy by 1.e4)
           qo3(i,k) = h1m4 * (do3v(i)+tthan(i)*(do3vp(i)-do3v(i)))
        END DO
     END DO
  END IF
  DO  i=1,ipts
     !visible and near ir diffuse albedo
     alvd(i) = sfcalb(i)
     alnd(i) = sfcalb(i)
     !visible and near ir direct beam albedo
     alvb(i) = sfcalb(i)
     alnb(i) = sfcalb(i)
     !visible and near ir direct beam albedo,if not ocean nor snow
     !function of cosine solar zenith angle..
     IF (slmsk(i+ibeg-1) == 0.0) THEN
        IF (sfcalb(i) <= 0.5) THEN
           albd0 = -18.0 * (0.5 - ACOS(coszen(i))/pi)
           albd0 = EXP (albd0)
           alvd1 = (alvd(i) - 0.054313) / 0.945687
           alnd1 = (alnd(i) - 0.054313) / 0.945687
           alvb(i) = alvd1 + (1.0 - alvd1) * albd0
           alnb(i) = alnd1 + (1.0 - alnd1) * albd0
        END IF
     END IF
  END DO
  !surface values of rrcl and ttcl
  DO  n=1,2
     DO  i=1,imax
        rrcl(i,n,1)=alvd(i)
        ttcl(i,n,1)=zero
     END DO
  END DO
  DO  n=3,nb
     DO  i=1,imax
        rrcl(i,n,1)=alnd(i)
        ttcl(i,n,1)=zero
     END DO
  END DO
  !end of cloud section  

  !the following code converts rrvco2,the volume mixing ratio of co2
  !into rrco2,the mass mixing ratio.
  rrvco2=rco2
  rrco2=rrvco2*ratco2mw
  IF(itimlw /= 0) THEN
     !long wave radiation 
     !account for reduced emissivity of any clds
     DO  k=1,lp1
        DO  i=1,imax
           eqcmt(i,k)=camt(i,k)*emcld(i,k)
        END DO
     END DO
     !get cld factor for lw calculations
     CALL clo89(cldfac,eqcmt,nclds,kbtm,ktop)
     !long wave radiation
     CALL lwr88(heatra,grnflx,topflx, press,temp,rh2o,qo3,cldfac,  &
          eqcmt,nclds,ktop,kbtm)
     DO  i=1,ipts
        ir = i + ibeg - 1
        flwup(ir) = topflx(i) * .001e0
        grnflx(i)=q14330*(hsigma*temp(i,lp1)**4-grnflx(i))
        !get lw flux down and up at ground(watts/m**2) - grnflx=lw down.
        flwdns(ir)=grnflx(i)/(1.43306e-06*1000.e0)
        flwups(ir)=hsigma*.001e0 * temp(i,lp1)**4
     END DO
     !convert heating rates to deg/sec
     DO  k=1,l
        DO  i=1,ipts
           hlw(i+ibeg-1,k)=heatra(i,k)*daysec
        END DO
     END DO
  END IF
  IF(itimsw /= 0) THEN
     !sw
     CALL swr93(fsw,hsw,uf,df,fswl,hswl,ufl,dfl,press,coszen,tauda,rh2o,&
          rrco2,ssolar,qo3, nclds,ktop,kbtm,camt,rrcl,ttcl,alvb, &
          alnb,alvd,alnd,gdfvb,gdfnb,gdfvd,gdfnd)

     !get sw fluxes in watts/m**2
     DO  i=1,ipts
        ir = i + ibeg - 1
        fswup(ir) = uf(i,1) * 1.e-3
        fswdn(ir) = df(i,1) * 1.e-3
        fswups(ir) = uf(i,lp1) * 1.e-3
        !c..couple w/m2 diff, if fswdns(ir)=df(i,lp1)*1.#e-3
        fswdns(ir) = (gdfvb(i)+gdfnb(i)+gdfvd(i)+gdfnd(i)) * 1.e-3
        !downward sfc flux for the sib parameteration
        !visible and near ir diffuse
        gdfvdr(ir) = gdfvd(i) * 1.e-3
        gdfndr(ir) = gdfnd(i) * 1.e-3
        !visible and near ir direct beam
        gdfvbr(ir) = gdfvb(i) * 1.e-3
        gdfnbr(ir) = gdfnb(i) * 1.e-3
     END DO
     !convert heating rates to deg/sec
     DO  k=1,l
        DO  i=1,ipts
           swh(i+ibeg-1,k)=hsw(i,k)*daysec
        END DO
     END DO
  END IF
1000 FORMAT(1H ,' YOU ARE CALLING GFDL RADIATION CODE FOR',i5,' PTS',  &
       'AND',i4,' LYRS,WITH KDAPRX,KO3,KCZ,KEMIS,KALB = ',5I2)
END SUBROUTINE radfs
