SUBROUTINE radtn
  !
  ! $Author: CPTEC $
  ! $Date: 2000/02/24 13:15:43 $
  ! $Revision: 1.0 $
  !
  !     ******************************************************************
  !$$$  subprogram documentation block
  !                .      .    .
  ! subprogram:    radtn       the outer radiation driver
  !   prgrmmr: black           org: w/np2      date: 93-12-??
  ! abstract:
  !     radtn primarily serves to set up the arrays needed as input
  !     for radfs (the inner radiation driver).  groups of model columns
  !     are sent to radfs but first they are "lifted" so that the lowest
  !     layer above the ground has a vertical index value of lm not lmh.
  !     this routine is called as often as desired (every 1 to 2 hours)
  !     for both the short and longwave effects.  the resulting temper-
  !     ature tendencies, total downward and shortwave upward fluxes are
  !     collected.
  !     the initial ground potential temperature is also computed here
  !     and is simply an adiabatic extrapolation from the lowest mid-
  !     layer value above the ground.
  ! program history log:
  !   87-09-??  black      - originator
  !   92-10-??  baldwin    - various cloud effects were included
  !                          which were already in the mrf
  !   93-11-??  zhao       - tied to updated gfdl radiation scheme
  !                          using model-predicted cloud
  !   95-03-25  black      - conversion from 1-d to 2-d in horizontal
  !   95-04-13  black      - parallelized the large loop stepping
  !                          through the domain that calls radfs
  !   95-10-10  zhao       - i) the calculation of cloud fraction was
  !                             changed to use both cloud water/ice
  !                             mixing ratio and relative humidity
  !                             (randall, 1994);
  !                          ii) the cloud inputs were changed to use
  !                              cloud fraction in each model layer
  !                              after y.t. hou (1995).
  !   96-06-03  zhao       - snow albedo is changed according to
  !                          suggestions from ken mitchell and fei chen
  !   96-07-23  zhao       - add call to solard to calculate the non-
  !                          dimensional sun-earth distance r1 which
  !                          will be used in radfs to compute solar
  !                          constant solc on each day
  !   96-07-26  black      - added ozone computations
  !
  ! usage: call radtn from main program ebu
  !
  ! note: convective clouds are added in this subroutine         
  ! for use in the eta model in which model-predicted      
  ! clouds are not used in the convective precipitation    
  ! processes.                                             
  ! for use with the version of the eta model in which     
  ! the model-predicted clouds are linked into the model's 
  ! convective precipitation processes, just set:          
  ! cncld=.false.                               
  !                               qingyun  zhao  12-9-94     

  IMPLICIT NONE

  INCLUDE "parmeta.f90"
  INCLUDE "parm.tbl.f90"

  REAL, PARAMETER     :: lda=lm+9,la=13
  REAL, PARAMETER     :: h1=1.e0,h2=2.e0,h5=5.e0,h10=10.e0,h1e5=1.e5,h1em5=1.e-5
  REAL, PARAMETER     :: d00=.00e0,d125=.125e0,d25=.25e0,d50=.50e0,d608=.608e0
  REAL, PARAMETER     :: d80=.80e0,capa=0.28589641e0,rtd=57.2957795e0 , wa=.10e0
  REAL, PARAMETER     :: wg=h1-wa,ksmud=0
  REAL, PARAMETER     :: a1=610.78e0,a2=17.2693882e0,a3=273.16e0,a4=35.86e0
  REAL, PARAMETER     :: pq0=379.90516e0,snoalb=0.55
  REAL, PARAMETER     :: slpm=1.01325e5,epsq1=1.e-5,epsq=2.e-12,epso3=1.e-10
  REAL, PARAMETER     :: hpinc=1.e1,cldrh0=0.80,tresh=1.00,rnrm=1./(tresh-cldrh0)
  REAL, PARAMETER     :: cldrh2=0.90,tresh2=1.00,rnrm2=1./(tresh2-cldrh2) 
  REAL, PARAMETER     :: clapse=-0.0005,clpse=-0.0006,dclps=-0.0001  
  REAL, PARAMETER     :: cm1=2937.4,cm2=4.9283,cm3=23.5518,eps=0.622,pbot=10000.0
  REAL, PARAMETER     :: stbol=5.67e-8,pi2=2.*3.14159265,rlag=14.8125
  INTEGER, PARAMETER :: nb=12,nsoil=2,nroot=2
  INTEGER, PARAMETER :: imjm=im*jm-jm/2,imt=2*im-1,jam=6+2*(jm-10)
  INTEGER, PARAMETER :: lm1=lm-1,lp1=lm+1

  LOGICAL               :: call1,short,long
  LOGICAL               :: bitx,bity,bitz,bitw,bit1,bit2,bitc
  LOGICAL               :: bits,bitcp1,bitsp1,cncld
  REAL    :: cfo3
  REAL, DIMENSION(3)  :: cthk, taucv
  REAL, DIMENSION(4)  :: ptopc
  REAL, DIMENSION(nb) :: abcff, pwts
  INTEGER, DIMENSION(im,jm) :: lvl
  INTEGER,DIMENSION(3) :: ltop

  INCLUDE "ctlblk.f90"
  INCLUDE "loops.f90"
  INCLUDE "masks.f90"
  INCLUDE "dynamd.f90"
  INCLUDE "phys.f90"
  INCLUDE "vrbls.f90"
  INCLUDE "pvrbls.f90"
  INCLUDE "soil.f90"
  INCLUDE "cldwtr.f90"
  INCLUDE "cnvcld.f90"
  INCLUDE "indx.f90"
  INCLUDE "acmcld.f90"
  INCLUDE "acmrdl.f90"
  INCLUDE "acmrds.f90"

  COMMON /swrsav/abcff,pwts,cfco2,cfo3,reflo3,rrayav
  COMMON /rd1tim/k400,cthk,ltop,ptopc,taucv,lvl

  REAL, DIMENSION(lm) :: tendk,cldamt
  REAL, DIMENSION(imt):: albedo,xlat,cosz,psfc,tskn,slmsk,umid,bmid
  REAL, DIMENSION(imt):: cv,flwup,fswdn,fswup,cemis,cstr
  REAL, DIMENSION(imt):: fswdns,fswups,flwdns,flwups,tauc,cvb,cvt,taudar
  REAL, DIMENSION(imt,3):: cldfrc
  REAL, DIMENSION(imt,lp1):: cldf,camt,pint,emis
  REAL, DIMENSION(imt,lm) :: tendl,pmid,tmid,tends,qmid,thmid,ozn,pozn
  REAL, DIMENSION(imt,lm) :: ccr,csmid,wmid,hmid,ccmid
  REAL, DIMENSION(im,jm)  ::  pdsl,pbotl,ptopl,pbotm,ptopm,pboth,ptoph,tot
  REAL, DIMENSION(0:im+1,0:jm+1) :: fne,fse,tl
  REAL, DIMENSION(9)  :: cc,ppt
  REAL, DIMENSION(lp1) :: phalf
  REAL, DIMENSION(imt,nb,lp1) :: rrcl,ttcl  
  INTEGER, DIMENSION(imt,lp1) :: ityp,ktop,kbtm
  INTEGER, DIMENSION(imt) :: icvb,icvt,nclds,kcld,kcvb,kcvt,numtyp
  INTEGER, DIMENSION(imt,3):: mtop,mbot
  INTEGER, DIMENSION(imt,lm) :: iw
  LOGICAL, DIMENSION(imt) :: bcld,btemp1

  INTEGER :: i,j,l,n,lml,lin,ivi,ir,ii,jd,ll,nfile,k400,kcco2,ktop1,kbt2,kbt1
  INTEGER :: nradpp,ntsph,itimsw,itimlw,iwkl,l400,ltrop,lbase,nc,nmod,kth2,kth1
  INTEGER :: lltop,llbot,malvl,nlvl,lmin1,nband,nbtm,nktp,ncld,kl,ks
  REAL :: plomd,pmdhi,phitp,p400,plbtm,ctim,ctim1,ctim2,utim,cfco2,cfco3,reflo3
  REAL :: rrayav,pm1,apem,tsfck,psfck,psum,pdif,pss,r1,fiw,fiq,u00kl,qint,qi
  REAL :: apem1, time, dayi, addl,rang,rcos1,rcos2,rsin1,times,exner,qw,pp,bi,ai 
  REAL :: pmod,clpfil,dthdp,ddp,rqkl,qc,tmt15,tmt0,cwmkl,qkl,tkl,hh,climit
  REAL :: cclimit,us,snofac,apes,clstp,dpc,cr1,cl2,cl1,sumd,cldmeb,cldmax,dd
  REAL :: tcld,delp,prs2,prs1,qsum,p2,p1,cc2,cc1,dpcl,hour,aa,ee,ff,bb,gg,denom
  REAL :: fctra,fctrb,pdslij,cfravg

  DATA plomd/64200.e0/,pmdhi/35000.e0/,phitp/15000.e0/,p400/40000.e0/  &
       , ltop/0,0,0/,plbtm/105000.e0/
  DATA nfile/14/
  DATA cc/0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8/
  DATA ppt/.14,.31,.70,1.6,3.4,7.7,17.,38.,85./
  DATA cthk/20000.0,20000.0,20000.0/
  DATA taucv/0.16, 0.14, 0.12/  &
       , cemis/imt*1./
  SAVE r1 !save r1 for uses in each radiation call

  ctim1=0.
  ctim2=24.*3600.
  ctim =ntsd*dt
  IF(ctim < ctim1)THEN
     utim=0.
  ELSE
     IF(ctim <= ctim2)THEN
        utim=(ctim-ctim1)/(ctim2-ctim1)
     ELSE
        utim=1.
     END IF
  END IF
  !convection
  cncld=.TRUE.
  !assign the pressures for cloud domain boundaries
  ptopc(1)=plbtm
  ptopc(2)=plomd
  ptopc(3)=pmdhi
  ptopc(4)=phitp

  !find the 'sea level pressure'.
  DO  j=1,jm
     DO  i=1,im
        pdsl(i,j)=res(i,j)*pd(i,j)
     END DO
  END DO

  !carry out the one-time calculations.
  !if(ntsd > 1)go to 110 !alterado por luiz flavio 16/12/99
  IF(ntsd<=1) THEN
     IF(.NOT.restrt) THEN
        !at this point, we must calculate the initial potential temperature
        !of the surface and of the subground.
        !extrapolate down for initial surface potential temperature.
        !mic$ do all autoscope
        DO  j=1,jm
           DO  i=1,im
              pm1=pdsl(i,j)*aeta(lmh(i,j))+pt
              apem1=(h1e5/pm1)**capa
              ths(i,j)=t(i,j,lmh(i,j))*(h1+d608*q(i,j,lmh(i,j)))*apem1
              tsfck=t(i,j,lmh(i,j))*(h1+d608*q(i,j,lmh(i,j)))
              psfck=pd(i,j)+pt
              IF(sm(i,j) < d50) THEN
                 qs(i,j)=pq0/psfck*EXP(a2*(tsfck-a3)/(tsfck-a4))
              ELSE IF(sm(i,j) > d00) THEN
                 ths(i,j)=sst(i,j)*(h1e5/(pd(i,j)+pt))**capa
              END IF
              tg(i,j)=(ths(i,j)*wa+tg(i,j)*wg)*(h1-sice(i,j)) +tg(i,j)*sice(i,j)
           END DO
        END DO
     END IF

     !calculate the number of steps at each point.
     !the array 'lvl' will coordinate vertical locations between
     !the lifted working arrays and the fundamental model arrays.
     !lvl holds the height (in model layers) of the topography at
     !each grid point.
     DO  j=1,jm
        DO  i=1,im
           lvl(i,j)=lm-lmh(i,j)
        END DO
     END DO

     !determine model layer limits for high(3), middle(2),
     !and low(1) clouds.  also find model layer that is just below
     !(height-wise) 400 mb. (k400)
     k400=0
     psum=pt
     pdif=slpm-pt
     DO  l=1,lm
        psum=psum+deta(l)*pdif
        IF(ltop(3) == 0)THEN
           IF(psum > phitp)ltop(3)=l
        ELSE IF(ltop(2) == 0)THEN
           IF(psum > pmdhi)ltop(2)=l
        ELSE IF(k400 == 0)THEN
           IF(psum > p400)k400=l
        ELSE IF(ltop(1) == 0)THEN
           IF(psum > plomd)THEN
              ltop(1)=l
              EXIT
           END IF
        END IF
     END DO

     !   105 continue !comentado por luiz flavio em 16/12/99

     !call gradfs once to calc. constants and get o3 data
     kcco2=0

     !calculate the midlayer pressures in the standard atmosphere
     pss=101325.
     pdif=pss-pt
     DO  l=1,lm1
        phalf(l+1)=aeta(l)*pdif+pt
     END DO
     phalf(1)=0.
     phalf(lp1)=pss

     CALL gradfs(phalf,kcco2,nfile)

     !call solard to calculate non-dimensional sun-earth distance
     CALL solard(r1)

     !call zenith simply to get the day of the year for
     !the setup of the ozone data
     time=(ntsd-1)*dt
     CALL zenith(time,dayi,hour)
     addl=0.
     IF(MOD(idat(3),4) == 0)addl=1.
     rang=pi2*(dayi-rlag)/(365.25+addl)
     rsin1=SIN(rang)
     rcos1=COS(rang)
     rcos2=COS(2.*rang)
     CALL o3clim
  END IF
  !110 continue

  !the following code is executed each time the radiation is called.
  !--convection--
  !  nradpp is the number of time steps to accumulate convective precip
  !  for radiation.
  !  note: this will not work if nrads and nradl are different unless
  !         they are integer multiples of each other
  !  clstp is the number of hours of the accumulation period
  ntsph=NINT(3600./dt)
  nradpp=MIN(nrads,nradl)
  clstp=1.0*nradpp/ntsph

  !convection
  !state whether the short or longwave computations are to be done.
  short=.FALSE.
  long=.FALSE.
  IF(MOD(ntsd,nrads) == 1.OR.restrt)short=.TRUE.
  IF(MOD(ntsd,nradl) == 1.OR.restrt)long=.TRUE.
  itimsw=0
  itimlw=0
  IF(short)itimsw=1
  IF(long) itimlw=1

  !find the mean cosine of the solar zenith angle
  !between the current time and the next time radiation is
  !called.  only average if the sun is above the horizon.
  time=(ntsd-1)*dt
  CALL zenith(time,dayi,hour)
  jd=INT(dayi+d50)
  IF(short)THEN
     !mic$ do all autoscope
     DO j=1,jm
        DO i=1,im
           czmean(i,j)=0.
           tot(i,j)=0.
        END DO
     END DO
     DO ii=0,nrads,nphs
        times=(ntsd-1)*dt+ii*dt
        CALL zenith(times,dayi,hour)
        !mic$ do all autoscope
        DO j=1,jm
           DO i=1,im
              IF(czen(i,j) > 0.)THEN
                 czmean(i,j)=czmean(i,j)+czen(i,j)
                 tot(i,j)=tot(i,j)+1.
              END IF
           END DO
        END DO
     END DO
     !mic$ do all autoscope
     DO j=1,jm
        DO i=1,im
           IF(tot(i,j) > 0.)czmean(i,j)=czmean(i,j)/tot(i,j)
        END DO
     END DO
  END IF

  !2345678901234567890123456789012345678901234567890123456789012345678901
  !mic$ parallel shared (acfrcv,acfrst,aeta,alb,capa,cc,cfrach)
  !mic$*        shared (cfracl,cfracm,clapse,clpse,clstp,cncld,cthk)
  !mic$*        shared (cuppt,czmean,d00,d50,dclps,epso3,epsq)
  !mic$*        shared (etad,glat,h1,h1e5,h1em5,h2,hbot,hour,hpinc,htop)
  !mic$*        shared (im,imjm,imt,itimlw,itimsw,jd,jm,k400,lm,lm1,lmh)
  !mic$*        shared (long,ltop,lvl,ncfrcv,ncfrst,pbot,pboth,pbotl)
  !mic$*        shared (pbotm,pd,pdsl,ppt,pt,ptopc,ptoph,ptopl,ptopm,q)
  !mic$*        shared (rcos1,rcos2,rlwin,rlwout,rlwtoa,rlwtt,rsin1)
  !mic$*        shared (rswin,rswout,rswtoa,rswtt,rtd,short,sigt4,sm)
  !mic$*        shared (sno,snoalb,stbol,t,taucv,ths)
  !mic$*        shared (hbm2,u00,cwm,htm,pq0,a2,a3,a4,utim,ul)
  !mic$*        shared (epsq1,lp1,nb,abcff)
  !mic$*        shared (vegfrc,irad,r1)
  !mic$*       private (albedo,apes,bcld,bit1,bit2,bitw,bitx,bity,bitz)
  !mic$*       private (btemp1,cc1,cc2,cemis,cfravg,cl1,cl2,cldamt)
  !mic$*       private (cldfrc,cldmax,cldmeb,clpfil,clstr,cmid,cosz,cr1)
  !mic$*       private (cstr,cv,delp,denom,dpcl,dthdp,emis,exner,fctra)
  !mic$*       private (fctrb,flwdns,flwup,flwups,fswdn,fswdns,fswup)
  !mic$*       private (fswups,i,icvb,icvt,ir,ivi,j,kbot,kbt1,kbt2,kcvb)
  !mic$*       private (kcvt,kl,kmid,kth1,kth2,ktop,ktop1,l,lbase,lin,ll)
  !mic$*       private (llbot,lltop,lmin,lmin1,lml,ltrop,malvl,mbot,mtop)
  !mic$*       private (n,nc,nlvl,nmod,numtyp,ozn,p1,p2,pdslij,pint)
  !mic$*       private (pmid,pmod,pozn,psfc,qmid,slmsk,sumd,tau0,tauc)
  !mic$*       private (taudar,tcld,tendl,tends,thmid,tmid,tskn,xlat)
  !mic$*       private (csmid,wmid,ccmid,iw,ccr,hmid,bmid,umid,us)
  !mic$*       private (cclimit,climit,hh,tkl,qkl,cwmkl,tmt0,tmt15)
  !mic$*       private (ai,bi,pp,qw,qi,qint,u00kl,fiq,iwkl,fiw,qc)
  !mic$*       private (rqkl,ddp,l400,kbtm,camt,ityp,kcld,nband,rrcl,ttcl)
  !mic$*       private (bitc,bits,bitcp1,bitsp1,nclds,ncld)
  !mic$*       private (qsum,nktp,nbtm,prs1,prs2,dd,ee,ff,aa,bb,gg)
  !mic$*       private (snofac,ir,pdslij)

  !********************************************************************
  !***  statement 120 is the beginning of the j loop through the domain
  !********************************************************************
  !mic$ do parallel
  DO  j = 1, jm-2, 2
     DO  l=1,lm
        DO  i=1,imt
           ir=irad(i)
           tmid(i,l)=t(ir,j,1)
           qmid(i,l)=epsq
           csmid(i,l)=d00
           wmid(i,l)=d00
           ccmid(i,l)=d00
           iw(i,l)=d00
           ccr(i,l)=d00
           hmid(i,l)=d00
           ozn(i,l)=epso3
           tends(i,l)=d00
           tendl(i,l)=d00
        END DO
     END DO
     DO n=1,3
        DO  i=1,imt
           cldfrc(i,n)=d00
           mtop(i,n)=0
           mbot(i,n)=0
        END DO
     END DO

     !fill in working arrays where values at l=lm are those that
     !are actually at eta level l=lmh.
     DO  i=1,imt
        ir=irad(i)
        lml=lmh(i,j)
        lin=lvl(i,j)
        bmid(i)=hbm2(i,j)
        umid(i)=u00(i,j)
        DO  l=1,lml
           pmid(i,l+lin)=aeta(l)*pdsl(i,j)+pt
           pint(i,l+lin+1)=etad(l+1)*pdsl(i,j)+pt
           exner=(h1e5/pmid(i,l+lin))**capa
           tmid(i,l+lin)=t(ir,j,l)
           thmid(i,l+lin)=t(ir,j,l)*exner
           qmid(i,l+lin)=q(ir,j,l)
           wmid(i,l+lin)=cwm(i,j,l)
           hmid(i,l+lin)=htm(i,j,l)
        END DO

        !fill in artificial values above the top of the domain.
        !pressure depths of these layers is 1 hpa.
        !temperatures above are already isothermal with (true) layer 1.
        IF(lin > 0)THEN
           DO  ivi=1,lin
              l=lin+1-ivi
              pmid(i,l)=pt-REAL(ivi+ivi-1)*0.5*hpinc
              pint(i,l+1)=pmid(i,l)+0.5*hpinc
              exner=(h1e5/pmid(i,l))**capa
              thmid(i,l)=tmid(i,l)*exner
           END DO
        END IF
        IF(lin == 0) THEN
           pint(i,1)=pt
        ELSE
           pint(i,1)=pmid(i,1)-0.5*hpinc
        END IF
     END DO

     !fill in the surface pressure, skin temperature, geodetic latitude,
     !zenith angle, sea mask, and albedo.  the skin temperature is
     !negative over water.
     DO  i=1,imt
        psfc(i)=pd(i,j)+pt
        apes=(psfc(i)*h1em5)**capa
        tskn(i)=ths(i,j)*apes*(h1-h2*sm(i,j))
        slmsk(i)=sm(i,j)
        !the following part is changed by q. zhao according to the
        !suggestions from ken. mitchell and fei chen.
        !albedo modified by snow, k. mitchell 03/96
        sno(i,j)=amax1(sno(i,j),0.)
        snofac=amin1(sno(i,j)/0.02, 1.0)
        albedo(i)=alb(i,j)+(1.0-vegfrc(i,j))*snofac*(snoalb-alb(i,j))
        xlat(i)=glat(i,j)*rtd
        cosz(i)=czmean(i,j)
     END DO

     !stratiform cloud section
     !--cloud--
     !  calculate stratiform cloud coverage at each model grid point which
     !  will be used in the model radiation           parameterization scheme.
     us=h1
     cclimit=1.0e-3
     climit =1.0e-20

     !qw, qi and qint
     DO  i=1,imt
        lml=lmh(i,j)
        lin=lvl(i,j)
        DO  l = 1, lml
           ll=l+lin
           hh=hmid(i,ll)*bmid(i)
           tkl=tmid(i,ll)
           qkl=qmid(i,ll)
           cwmkl=wmid(i,ll)
           tmt0=(tkl-273.16)*hh
           tmt15=amin1(tmt0,-15.)*hh
           ai=0.008855
           bi=1.
           IF(tmt0 < -20.)THEN
              ai=0.007225
              bi=0.9674
           END IF
           pp=pmid(i,ll)
           qw =hh*pq0/pp*EXP(hh*a2*(tkl-a3)/(tkl-a4))
           qi=qw *(bi+ai*amin1(tmt0,0.))
           qint=qw *(1.-0.00032*tmt15*(tmt15+15.))
           IF(tmt0 <= -40.) qint=qi
           !ice-water id number iw
           u00kl=umid(i)+ul(l)*(0.95-umid(i))*utim
           IF(tmt0 < -15.0) THEN
              fiq=qkl-u00kl*qi
              IF(fiq > d00.OR.cwmkl > climit) THEN
                 iw(i,ll)=1
              ELSE
                 iw(i,ll)=0
              END IF
           END IF

           IF(tmt0 >= 0.0) THEN
              iw(i,ll)=0
           END IF
           IF(tmt0 < 0.0.AND.tmt0 >= -15.0) THEN
              iw(i,ll)=0
              IF(iw(i,ll-1) == 1.AND.cwmkl > climit) iw(i,ll)=1
           END IF
           iwkl=iw(i,ll)
           !the satuation specific humidity------------------------
           fiw=float(iwkl)
           qc    =(h1-fiw)*qint+fiw*qi
           !the relative humidity----------------------------------
           IF(qc <= epsq1.OR.qkl <= epsq1) THEN
              rqkl=d00
           ELSE
              rqkl=qkl/qc
           END IF
           !cloud cover ratio ccr----------------------------------
           IF(rqkl >= 0.9999)       THEN
              ccr(i,ll)=amin1(us,rqkl)
           ELSE
              ccr(i,ll)= rqkl*(h1-EXP((-1000.0*cwmkl)/(us-rqkl)))
           END IF
           csmid(i,ll)=ccr(i,ll)
        END DO
     END DO
     !--cloud--
     ! now check the clouds produced above to make sure they are good
     ! enough for radiation calculations

     !   no stratiform clouds for this type
     DO  i=1,imt
        lml=lmh(i,j)
        lin=lvl(i,j)
        !zero out cldamt if land and below pbot above ground
        IF(sm(i,j) < 0.5) THEN
           DO   l=1,lml
              ll=lml-l+1+lin
              ddp=psfc(i)-pmid(i,ll)
              IF(ddp >= pbot) EXIT
              csmid(i,ll)=d00
           END DO
        END IF

        !check for ocean stratus (low cloud)
        !look only over ocean and only if an inversion (dthdp.le.-0.05)
        !is present with at least 2 cloud free layers above it
        IF(sm(i,j) > 0.5) THEN
           !find base of inversion
           lbase=lm
           DO  l=1, lml-1
              ll=lml-l+1+lin
              dthdp=(thmid(i,ll-1)-thmid(i,ll)) /(pmid(i,ll-1)-pmid(i,ll))
              IF (dthdp <= clapse) THEN
                 lbase=ll
                 EXIT
              END IF
           END DO
           IF(csmid(i,lbase-1) <= d00.AND.csmid(i,lbase-2) <= d00  &
                .AND.lbase < lm) THEN
              IF (dthdp > clpse) THEN
                 clpfil = 1.0 - ((clpse - dthdp) / dclps)
              ELSE
                 clpfil = 1.0
              END IF
              DO  l=1,lml
                 ll=lml-l+1+lin
                 ddp=psfc(i)-pmid(i,ll)
                 IF(ddp >= pbot) EXIT
                 csmid(i,ll)=csmid(i,ll)*clpfil
              END DO

              !  if no inversion or if clds exist in either of the 2 layers above
              !  inversion, zero out cloud below lmin1 (pbot)
           ELSE
              DO   l=1,lml
                 ll=lml-l+1+lin
                 ddp=psfc(i)-pmid(i,ll)
                 IF(ddp >= pbot) EXIT
                 csmid(i,ll)=d00
              END DO
           END IF
        END IF
        !remove high clouds above the tropopause
        l400=lm
        DO  l = 1, lml
           ll=lml-l+1+lin
           IF(pmid(i,ll) <= 40000.0) THEN
              l400=ll
              EXIT
           END IF
        END DO
        ltrop=lm
        DO  ll=l400,2,-1
           dthdp=(thmid(i,ll-1)-thmid(i,ll)) /(pmid(i,ll-1)-pmid(i,ll))
           IF (dthdp < -0.0025.OR.qmid(i,ll) <= epsq1) THEN
              ltrop=ll
              EXIT
           END IF
        END DO
        IF (ltrop < lm) THEN
           DO  ll=ltrop,1,-1
              csmid(i,ll)=d00
           END DO
        END IF
     END DO

     !determine the fractional cloud coverage for high, mid
     !and low of clouds from the cloud coverage at each level
     !note: this is for diagnostic only!!!

     !
     ! essa area eh prioritaria. deve-se forca-la a nao vetorizacao
     ! caso o problema persistir.
     !
     loop300: DO  i=1,imt ! 
        !--convection--                                             
        !if there will be convective cloud, skip the stratiform     
        !cloud computation                                          
        IF(cncld) THEN                      
           pmod=cuppt(i,j)*24.0*1000.0/clstp 
           IF (pmod > ppt(1)) CYCLE       !loop300 
        END IF
        cstr(i)=0.0
        !establish the lowest cloud bottom (lmin1) as the lowest
        !(mid-) level that is at least pbot above the surface
        DO  l=1,lm1
           ll=lm1-l+1
           lmin1=ll
           IF (psfc(i)-pmid(i,ll) >= pbot) EXIT                      
        END DO
        DO  nlvl=1,3                  
           cldmax=d00                   
           cldmeb=d00                 
           sumd  =d00                  
           malvl =lm                   
           lltop=ltop(nlvl)+lvl(i,j)    

           !go to the next point if the top of the cloud-type in question 
           !is below ground or is in the lowest layer above ground.       
           IF(lltop >= lm)CYCLE loop300  

           IF(nlvl > 1)THEN
              llbot=ltop(nlvl-1)-1+lvl(i,j)
              llbot=MIN(llbot,lm1)
           ELSE
              llbot=lm1
           END IF
           DO  l=lltop,llbot
              cldamt(l)=csmid(i,l)
              IF(cldamt(l) > cldmax) THEN
                 malvl=l
                 cldmax=cldamt(l)
              END IF
           END DO

           ! now, calculate the total cloud fraction in this pressure domain
           ! using the method developed by y.h., k.a.c. and a.k. (nov., 1992).
           ! in this method, it is assumed that seperated cloud layers are
           ! radomly overlapped and adjacent cloud layers are maximum overlapped.
           ! vertical location of each type of cloud is determined by the thickest
           ! continuing cloud layers in the domain.
           !*********************************************************************
           cl1=0.0
           cl2=0.0
           kbt1=llbot
           kbt2=llbot
           kth1=0
           kth2=0
           DO  ll=lltop, llbot
              l=llbot-ll+lltop
              bit1=.FALSE.
              cr1=cldamt(l)
              bitx=(pint(i,l) >= ptopc(nlvl+1)).AND.  &
                   (pint(i,l) < ptopc(nlvl)).AND. (cldamt(l) > 0.0)
              bit1=bit1.OR.bitx
              IF(.NOT.bit1) CYCLE
              !bity=t: first cloud layer; bitz=t:consecutive cloud layer
              !note:  we assume that the thickness of each cloud layer in the
              !       domain is less than 200 mb to avoid too much cooling or
              !       heating. so we set cthk(nlvl)=200*e2. but this limit may
              !       work well for convective clouds. modification may be
              !       needed in the future.
              bity=bitx.AND.(kth2 <= 0)
              bitz=bitx.AND.(kth2 > 0)
              IF(bity) THEN
                 kbt2=l
                 kth2=1
              END IF
              IF(bitz) THEN
                 ktop1=kbt2-kth2+1
                 dpcl = pmid(i,kbt2)-pmid(i,ktop1)
                 IF(dpcl < cthk(nlvl)) THEN
                    kth2=kth2+1
                 ELSE
                    kbt2=kbt2-1
                 END IF
              END IF
              IF(bitx) cl2=amax1(cl2,cr1)

              !at the domain boundary or separated cld layers, random overlap.
              !choose the thickest or the largest fraction amt as the cld
              !layer in that domain.
              bit2=.FALSE.
              bity=bitx.AND.(cldamt(l-1) <= 0.0.OR. pint(i,l-1) < ptopc(nlvl+1))
              bitz=bity.AND.cl1 > 0.0
              bitw=bity.AND.cl1 <= 0.0
              bit2=bit2.OR. bity
              IF(.NOT.bit2) CYCLE
              IF(bitz) THEN
                 kbt1=INT((cl1*kbt1+cl2*kbt2)/(cl1+cl2))
                 kth1=INT((cl1*kth1+cl2*kth2)/(cl1+cl2))+1
                 cl1 =cl1+cl2-cl1*cl2
              END IF
              IF(bitw) THEN
                 kbt1=kbt2
                 kth1=kth2
                 cl1 =cl2
              END IF
              IF(bity) THEN
                 kbt2=llbot
                 kth2=0
                 cl2 =0.0
              END IF
           END DO
           cldfrc(i,nlvl)=cl1
           mtop(i,nlvl)=MIN(kbt1,kbt1-kth1+1)
           mbot(i,nlvl)=kbt1
        END DO
     END DO loop300

     !--convection--
     !convective cloud section
     ! this part was modified to compute convective clouds at each
     ! model layer based on convective precipitation rates. currently,
     ! clouds are set to 0.75*cv(i) below 400mb
     ! and    0.90*cv(i) above 400mb to account for cirrus cap
     !                                       q.zhao   95-3-22

     !   compute the convective cloud cover for radiation
     IF(cncld)  THEN
        DO  i=1,imt
           pmod=cuppt(i,j)*24.0*1000.0/clstp
           nmod=0
           DO  nc=1,9
              IF(pmod > ppt(nc)) nmod=nc
           END DO
           !cloud tops and bottoms come from cucnvc
           !add lvl to be consistent with other working arrays
           IF (nmod == 0) THEN
              cv(i)=d00
              icvt(i)=0
              icvb(i)=0
           ELSE IF (nmod == 9) THEN
              cv(i)=cc(9)
              icvt(i)=INT(htop(i,j)+d50)+lvl(i,j)
              icvb(i)=INT(hbot(i,j)+d50)+lvl(i,j)
           ELSE
              cc1=cc(nmod)
              cc2=cc(nmod+1)
              p1=ppt(nmod)
              p2=ppt(nmod+1)
              cv(i)=cc1+(cc2-cc1)*(pmod-p1)/(p2-p1)
              icvt(i)=INT(htop(i,j)+d50)+lvl(i,j)
              icvb(i)=INT(hbot(i,j)+d50)+lvl(i,j)
           END IF
        END DO
        !make sure clouds are deep enough
        DO  i=1,imt
           bcld(i) = cv(i) > d00.AND. (icvb(i) - icvt(i)) >= 2
           btemp1(i) = bcld(i)
        END DO
        !compute convective cloud fraction
        DO  i=1, imt
           IF(bcld(i)) THEN
              lml=lmh(i,j)
              lin=lvl(i,j)
              DO  l = 1, lml
                 ll=l+lin
                 IF(ll > icvb(i).OR.ll < icvt(i)) THEN
                    ccmid(i,ll)=d00
                 ELSE
                    IF(pmid(i,ll) < 40000.0) THEN
                       ccmid(i,ll)=0.90*cv(i)
                    ELSE
                       ccmid(i,ll)=0.75*cv(i)
                    END IF
                 END IF
              END DO
           END IF
        END DO
        !the following parts are for diagnostic only!!!
        !get cloud type of top and bottom
        DO  kl=1,3
           DO  i=1,imt
              lltop=ltop(kl)+lvl(i,j)
              IF (btemp1(i).AND. icvb(i) >= lltop) THEN
                 kcvb(i)=kl
                 btemp1(i)=.FALSE.
              END IF
           END DO
        END DO
        DO  i=1,imt
           bcld(i) = bcld(i) .AND. .NOT.btemp1(i)
           kcvt(i) = 0
           btemp1(i) = bcld(i)
        END DO
        DO  kl=1,3
           DO  i=1,imt
              lltop=ltop(kl)+lvl(i,j)
              IF (btemp1(i).AND. icvt(i) >= lltop) THEN
                 kcvt(i) = kl
                 btemp1(i) = .FALSE.
              END IF
           END DO
        END DO
        DO  i=1,imt
           bcld(i) = bcld(i) .AND. .NOT.btemp1(i)
           numtyp(i) = kcvt(i) - kcvb(i)
        END DO
        !different assumptions depending on how deep the convection is
        !numtyp=0,1,2 for 1,2,3 cloud domains covered by convection
        DO  i=1,imt
           IF(bcld(i)) THEN
              IF (numtyp(i) > 1) mtop(i,kcvb(i)) = icvb(i) - 1
              IF (numtyp(i) == 1.AND. kcvb(i) <= 1) mtop(i,kcvb(i)) = icvb(i) - 1
           END IF
        END DO
        DO  i=1,imt
           IF (bcld(i)) THEN
              !cumulus model - l = .75 * cc to account for low fractocumulus
              !(slingo, 1987)  m = .25 * cc and is tall convective tower
              !                h = .25 * cc ,but cirrus cap if cv (ie pmod)
              !                     is large enough and convection above 400mb
              IF (numtyp(i) > 1) THEN
                 !if convection is deep enough to span 3 types
                 !then proceed as advertised above (for l,m,h) lo=2 lyr thk
                 mbot(i,kcvb(i))   = icvb(i)
                 cldfrc(i,kcvb(i)) = 0.75*cv(i)
                 mtop(i,kcvt(i))   = icvt(i)
                 mbot(i,kcvt(i))   = icvt(i)
                 cldfrc(i,kcvt(i)) = 0.25*cv(i)
                 IF (cv(i) >= 0.40.AND.icvt(i) <= k400)  &
                      cldfrc(i,kcvt(i)) = 2.0 * (cv(i) - 0.3)
                 mtop(i,kcvb(i)+1)   = mbot(i,kcvt(i)) + 1
                 mbot(i,kcvb(i)+1)   = mtop(i,kcvb(i)) - 1
                 cldfrc(i,kcvb(i)+1) = 0.25*cv(i)
              ELSE IF (numtyp(i) == 1) THEN
                 !convection only spans 2 types
                 IF (icvb(i) <= 1) THEN
                    !if l,m clds
                    mbot(i,kcvb(i))   = icvb(i)
                    cldfrc(i,kcvb(i)) = 0.75*cv(i)
                    mtop(i,kcvt(i))   = icvt(i)
                    mbot(i,kcvt(i))   = mtop(i,kcvb(i)) - 1
                    cldfrc(i,kcvt(i)) = 0.25*cv(i)
                    !even though only 2 cloud types, if deep enough convection, then
                    !use increased cloud to simulate some type of thick mid cloud
                    !(.25 - .65)
                    IF (cv(i) >= 0.40.AND.icvt(i) <= k400)  &
                         cldfrc(i,kcvt(i)) = cv(i) - 0.15
                 ELSE
                    !if m,h clds
                    mtop(i,kcvt(i))   = icvt(i)
                    mbot(i,kcvt(i))   = icvt(i)
                    cldfrc(i,kcvt(i)) = 0.25*cv(i)
                    !model layer containing 400 mb is saved in k400
                    IF (cv(i) >= 0.40.AND.icvt(i) <= k400)  &
                         cldfrc(i,kcvt(i)) = 2.0 * (cv(i) - 0.3)
                    mtop(i,kcvb(i))   = mbot(i,kcvt(i)) + 1
                    mbot(i,kcvb(i))   = icvb(i)
                    !even though only 2 cloud types, use cloud amount for lower type
                    !that is mean of thin tower and low fractocu
                    cldfrc(i,kcvb(i)) = 0.50*cv(i)
                 END IF
              ELSE IF (numtyp(i) == 0) THEN
                 !convection only spans one cloud type...kcvb(i)=kcvt(i)
                 IF (icvb(i) <= 1) THEN
                    !low cloud...assume fractocu as above
                    mtop(i,kcvb(i)) = icvt(i)
                    mbot(i,kcvb(i)) = icvb(i)
                    cldfrc(i,kcvb(i)) = 0.75 * cv(i)
                 ELSE IF (icvb(i) == 2) THEN
                    !mid cloud (not thin tower)
                    mtop(i,kcvb(i)) = icvt(i)
                    mbot(i,kcvb(i)) = icvb(i)
                    cldfrc(i,kcvb(i)) = 0.50 * cv(i)
                 ELSE IF (icvb(i) > 2) THEN
                    !high cloud (not anvil)
                    mtop(i,kcvb(i)) = icvt(i)
                    mbot(i,kcvb(i)) = icvb(i)
                    cldfrc(i,kcvb(i)) = 0.50 * cv(i)
                 END IF
              END IF
           END IF
        END DO

        !remove high clouds above the tropopause
        l400=lm
        DO  i = 1, imt
           lml=lmh(i,j)
           lin=lvl(i,j)
           DO  l = 1, lml
              ll=lml-l+1+lin
              IF(pmid(i,ll) <= 40000.0) THEN
                 l400=ll
                 EXIT
              END IF
           END DO
           ltrop=lm
           DO  ll=l400,2,-1
              dthdp=(thmid(i,ll-1)-thmid(i,ll)) /(pmid(i,ll-1)-pmid(i,ll))
              IF (dthdp < -0.0025.OR.qmid(i,ll) <= epsq1) THEN
                 ltrop=ll
                 EXIT
              END IF
           END DO
           IF (ltrop < lm) THEN
              DO  ll=ltrop,1,-1
                 ccmid(i,ll)=d00
              END DO
           END IF
        END DO
     END IF
     !end of convective cloud section

     !--convection--
     !set the un-needed taudar to one
     DO i = 1, imt
        taudar(i)=1.0
     END DO

     ! now, calculate the cloud radiative properties after davis (1982),
     ! harshvardhan et al (1987) and y.h., k.a.c. and a.k. (1993).
     ! update: the following parts are modified, after y.t.h. (1994), to
     !         calculate the radiative properties of clouds on each model
     !         layer. both convective and stratiform clouds are used
     !         in this calculations.
     !                                     qingyun zhao   95-3-22

     !initialize arrays for uses later
     DO  i = 1, imt
        lml=lmh(i,j)
        lin=lvl(i,j)
        !note: layer=1 is the surface, and layer=2 is the first cloud
        !      layer above the surface and so on.
        emis(i,1)=1.0
        ktop(i,1)=lp1
        kbtm(i,1)=lp1
        camt(i,1)=1.0
        ityp(i,1)=0
        kcld(i)=2
        DO  nband=1,nb
           rrcl(i,nband,1)=0.0
           ttcl(i,nband,1)=1.0
        END DO

        DO  l=2,lp1
           ityp(i,l)=0
           camt(i,l)=0.0
           ktop(i,l)=1
           kbtm(i,l)=1
           emis(i,l)=0.0
           DO  nband=1,nb
              rrcl(i,nband,l)=0.0
              ttcl(i,nband,l)=1.0
           END DO
        END DO

        !now calculate the amount, top, bottom and type of each cloud layer
        !cloud type=1: stratiform cloud
        !      type=2: convective cloud
        ! when both convective and stratiform clouds exist at the same point,
        ! select convective cloud (type=2),in other words, convective clouds
        ! have the higher priority than stratiform clouds.
        ! cloud layers are separated by:
        !      1. no-cloud layer
        !      2. different cloud type
        !note: there is only one convective cloud layer in one column.
        !ktop and kbtm are the top and bottom of each cloud layer in terms o
        !eta model level.
        DO  l = 2, lml
           ll = lml - l +1 + lin
           bitc=ccmid(i,ll) > 0.1
           bits=csmid(i,ll) > 0.1
           bitcp1=ccmid(i,ll+1) > 0.1
           bitsp1=csmid(i,ll+1) > 0.1
           bit1=bits.OR.bitc
           IF(bit1) THEN
              IF(ityp(i,kcld(i)) == 0) THEN
                 camt(i,kcld(i))=csmid(i,ll)
                 ityp(i,kcld(i))=1
                 kbtm(i,kcld(i))=ll
                 IF(bitc) THEN
                    camt(i,kcld(i))=ccmid(i,ll)
                    ityp(i,kcld(i))=2
                 END IF
              ELSE
                 IF(bitc) THEN
                    IF(bitcp1) THEN
                       camt(i,kcld(i))=amax1(camt(i,kcld(i)),ccmid(i,ll))
                    ELSE
                       kcld(i)=kcld(i)+1
                       camt(i,kcld(i))=ccmid(i,ll)
                       ityp(i,kcld(i))=2
                       ktop(i,kcld(i)-1)=ll+1
                       kbtm(i,kcld(i))=ll
                    END IF
                 ELSE
                    IF(bitcp1) THEN
                       kcld(i)=kcld(i)+1
                       camt(i,kcld(i))=csmid(i,ll)
                       ityp(i,kcld(i))=1
                       ktop(i,kcld(i)-1)=ll+1
                       kbtm(i,kcld(i))=ll
                    ELSE
                       camt(i,kcld(i))=amax1(camt(i,kcld(i)),csmid(i,ll))
                    END IF
                 END IF
              END IF
           ELSE
              IF(bitcp1.OR.bitsp1) THEN
                 kcld(i)=kcld(i)+1
                 ktop(i,kcld(i)-1)=ll+1
                 ityp(i,kcld(i))=0
                 camt(i,kcld(i))=0.0
              END IF
           END IF
        END DO
        !the real number of cloud layers is (the first is the groung;
        !the last is the sky):
        nclds(i)=kcld(i)-2
        ncld=nclds(i)

        !now calculate cloud radiative properties
        IF(ncld >= 1) THEN
           !note: the following calculations, the unit for pressure is mb!!!
           DO  nc = 2, ncld+1
              tauc(i)=0.0
              qsum=0.0
              nktp=lp1
              nbtm=0
              bitx=camt(i,nc) > 0.1
              nktp=MIN(nktp,ktop(i,nc))
              nbtm=MAX(nbtm,kbtm(i,nc))
              DO  ll = nktp, nbtm
                 IF(ll >= ktop(i,nc).AND.ll <= kbtm(i,nc).AND.bitx) THEN
                    prs1=pint(i,ll)*0.01
                    prs2=pint(i,ll+1)*0.01
                    delp=prs2-prs1
                    tcld=tmid(i,ll)-273.16
                    qsum=qsum+qmid(i,ll)*delp*(prs1+prs2) /(120.1612*SQRT(tmid(i,ll)))
                    !for convective cloud or startiform cloud with top above 500mb
                    IF (ityp(i,nc) == 2 .OR.pint(i,ktop(i,nc)) <= ptopc(3)) THEN
                       IF(tcld <= -10.0) THEN
                          tauc(i)=tauc(i)+delp*amax1(0.1e-3, 2.0e-6*(tcld+82.5)**2)
                       ELSE  
                          tauc(i)=tauc(i)+delp*amin1(0.08,6.949e-3*tcld+0.1)
                       END IF
                    ELSE
                       !for low and mid stratiform clouds
                       IF(tcld <= -20.0) THEN
                          tauc(i)=tauc(i)+delp*amax1(0.1e-3,2.56e-5* (tcld+82.5)**2)
                       ELSE
                          tauc(i)=tauc(i)+delp*0.1
                       END IF
                    END IF
                 END IF
              END DO
              IF(bitx) emis(i,nc)=1.0-EXP(-0.75*tauc(i))
              IF(qsum >= epsq1) THEN
                 DO  nband=1,nb
                    IF(bitx) THEN
                       dd=tauc(i)/(tauc(i)+abcff(nband)*qsum)
                       ee=1.0-dd
                       ff=1.0-dd*0.85
                       aa=amin1(50.0,SQRT(3.0*ee*ff)*tauc(i))
                       aa=EXP(-aa)
                       bb=ff/ee
                       gg=SQRT(bb)
                       dd=(gg+1.0)*(gg+1.0)-(gg-1.0)*(gg-1.0)*aa*aa
                       rrcl(i,nband,nc)=amax1(0.1e-5, (bb-1.0)*(1.0-aa*aa)/dd)
                       ttcl(i,nband,nc)=amax1(0.1e-5, 4.0*gg*aa/dd)
                    END IF
                 END DO
              END IF
           END DO
        END IF
     END DO

     ! compute ozone at midlayers 
     !modify pressures so that the entire column of ozone (to 0 mb)
     !is included in the model column even when pt > 0 mb
     DO l=1,lm
        DO i=1,imt
           denom=1./(pint(i,lp1)-pint(i,1))
           fctra=pint(i,lp1)*denom
           fctrb=-pint(i,1)*pint(i,lp1)*denom
           pozn(i,l)=pmid(i,l)*fctra+fctrb
        END DO
     END DO

     CALL ozon2d(imt,lm,pozn,xlat,rsin1,rcos1,rcos2,ozn)
     CALL radfs (imt,imt  &
          ,     psfc,pmid,pint,qmid,tmid,ozn,tskn,slmsk,albedo,xlat  &
          ,     camt,ityp,ktop,kbtm,nclds,emis,rrcl,ttcl ,     cosz,taudar,1,imt  &
          ,     1,0 ,     etad,aeta,itimsw,itimlw,jd,r1,hour,tends,tendl  &
          ,     flwup,fswup,fswdn,fswdns,fswups,flwdns,flwups)
     DO  i=1,imt
        pdslij=pdsl(i,j)
        pmod=cuppt(i,j)*24.0*1000.0/clstp
        cfracl(i,j)=cldfrc(i,1)
        cfracm(i,j)=cldfrc(i,2)
        cfrach(i,j)=cldfrc(i,3)

        !arrays acfrst and acfrcv accumulate average stratiform and
        !convective cloud fractions, respectively.  this information
        !is passed to the post processor via common block acmcld.
        cfravg   = amax1(cfracl(i,j),amax1(cfracm(i,j),cfrach(i,j)))
        IF(cncld) THEN
           IF (pmod <= ppt(1)) THEN
              acfrst(i,j) = acfrst(i,j) + cfravg
              ncfrst(i,j) = ncfrst(i,j) + 1
           ELSE
              acfrcv(i,j) = acfrcv(i,j) + cfravg
              ncfrcv(i,j) = ncfrcv(i,j) + 1
           END IF
        ELSE
           acfrst(i,j) = acfrst(i,j) + cfravg
           ncfrst(i,j) = ncfrst(i,j) + 1
        END IF
     END DO
     !collect atmospheric temperature tendencies due to radiation.
     !also collect the total sw and incoming lw radiation (w/m**2)
     !and convert to form needed for prediction of ths in surfce.
     DO  i=1,imt
        DO  l=1,lm
           ll=lvl(i,j)+l
           IF(short)rswtt(i,j,l)=tends(i,ll)
           IF(long) rlwtt(i,j,l)=tendl(i,ll)
           IF(ll == lm)CYCLE
        END DO
     END DO
     !sum the lw incoming and sw radiation (w/m**2) for radin.
     DO  i=1,imt
        IF(long)THEN
           sigt4(i,j)=stbol*tmid(i,lm)*tmid(i,lm)* tmid(i,lm)*tmid(i,lm)
        END IF
        !accumulate various lw and sw radiative fluxes for post
        !processor.  passed via common acmrdl and acmrds.
        IF(long)THEN
           rlwin(i,j)  = flwdns(i)
           rlwout(i,j) = flwups(i)
           rlwtoa(i,j) = flwup(i)
        END IF
        IF(short)THEN
           rswin(i,j)  = fswdns(i)
           rswout(i,j) = fswups(i)
           rswtoa(i,j) = fswup(i)
        END IF
     END DO
  END DO
  !calls to radiation this time step are complete.

  !mic$ end parallel

  !horizontal smoothing of temperature tendencies
  IF(short) THEN
     !mic$ do all autoscope
     !mic$* private (fne,fse,tl)
     DO  l=1,lm
        CALL zero2(tl)
        CALL zero2(fne)
        CALL zero2(fse)
        DO  ks=1,ksmud
           DO  j=1,jm
              DO  i=1,im
                 tl(i,j)=rswtt(i,j,l)*htm(i,j,l)
              END DO
           END DO
           DO  j=1,jm-1
              DO  i=1,im
                 fne(i,j)=(tl(i+ihe(j),j+1)-tl(i,j)) *htm(i,j,l)*htm(i+ihe(j),j+1,l)
              END DO
           END DO
           DO  j=2,jm
              DO  i=1,im
                 fse(i,j)=(tl(i+ihe(j),j-1)-tl(i,j)) *htm(i+ihe(j),j-1,l)*htm(i,j,l)
              END DO
           END DO
           DO  j=3,jm-2
              DO  i=1,im
                 tl(i,j)=(fne(i,j)-fne(i+ihw(j),j-1)+fse(i,j)-fse(i+ihw(j),j+1))  &
                      *hbm2(i,j)*d125+tl(i,j)
              END DO
           END DO
           DO  j=1,jm
              DO  i=1,im
                 rswtt(i,j,l)=tl(i,j)
              END DO
           END DO
        END DO
     END DO
  END IF

  IF(long)THEN
     !mic$ do all autoscope
     !mic$* private (fne,fse,tl)
     DO  l=1,lm
        CALL zero2(tl)
        CALL zero2(fne)
        CALL zero2(fse)
        DO  ks=1,ksmud
           DO  j=1,jm
              DO  i=1,im
                 tl(i,j)=rlwtt(i,j,l)*htm(i,j,l)
              END DO
           END DO
           DO  j=1,jm-1
              DO  i=1,im
                 fne(i,j)=(tl(i+ihe(j),j+1)-tl(i,j)) *htm(i,j,l)*htm(i+ihe(j),j+1,l)
              END DO
           END DO
           DO  j=2,jm
              DO  i=1,im
                 fse(i,j)=(tl(i+ihe(j),j-1)-tl(i,j)) *htm(i+ihe(j),j-1,l)*htm(i,j,l)
              END DO
           END DO
           DO  j=3,jm-2
              DO  i=1,im
                 tl(i,j)=(fne(i,j)-fne(i+ihw(j),j-1)+fse(i,j)-fse(i+ihw(j),j+1))  &
                      *hbm2(i,j)*d125+tl(i,j)
              END DO
           END DO
           DO  j=1,jm
              DO  i=1,im
                 rlwtt(i,j,l)=tl(i,j)
              END DO
           END DO
        END DO
     END DO
  END IF

  !reset cuppt,htop,hbot from this call to radtn
  IF(MOD(ntsd,nradpp) == 1)THEN
     DO  j=1,jm
        DO  i=1,im
           cuppt(i,j)=d00
           htop(i,j)=100.
           hbot(i,j)=d00
        END DO
     END DO
  END IF
END SUBROUTINE radtn
