!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    SUBROUTINE RADTN
!     ******************************************************************
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!                .      .    .
! SUBPROGRAM:    RADTN       THE OUTER RADIATION DRIVER
!   PRGRMMR: BLACK           ORG: W/NP22     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
!   97-05-19  ZHAO       - DIAGNOSTIC CLOUDS (LOW, MIDDLE, AND HIGH)
!                          ARE MODIFIED TO USE THE MAXIMUM OF
!                          CONVECTIVE AND STRATIFORM. THIS WILL REPLACE
!                          THE PREVIOUS SCHEME WHICH USES ONLY CONVECTIVE
!                          CLOUDS AT CONVECTIVE POINTS. THIS WILL
!                          AFFECT CFRACL, CFRACM, CFRACH, AND WILL
!                          AFFECT THE TOTAL CLOUD FRACTION CALCULATION
!                          IN THE POST PROCESSORS.
!   98-??-??  TUCCILLO   - ADDED PARALLELISM FOR CLASS VIII
!   98-10-27  BLACK      - PARALLELISM INTO NEWEST VERSION



! USAGE: CALL RADTN FROM MAIN PROGRAM EBU
!   INPUT ARGUMENT LIST:
!     NONE

!   OUTPUT ARGUMENT LIST:
!     NONE

!   OUTPUT FILES:
!     NONE

!   SUBPROGRAMS CALLED:

!     UNIQUE:
!        ZENITH
!        RADFS
!        GRADFS
!        SOLARD
!        O3CLIM
!        OZON2D

!     LIBRARY:
!        NONE

!   COMMON BLOCKS: CTLBLK
!                  LOOPS
!                  MASKS
!                  DYNAMD
!                  PHYS
!                  VRBLS
!                  PVRBLS
!                  CLDWTR
!                  CNVCLD
!                  CUINIT
!                  INDX
!                  ACMCLD
!                  ACMRDS
!                  ACMRDL

! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!   MACHINE : IBM SP
!$$$
!     ******************************************************************
!     *   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     *
!     ******************************************************************
!-----------------------------------------------------------------------
    INCLUDE "parmeta.f90"
    INCLUDE "parm.tbl.f90"
    INCLUDE "parmsoil.f90"
    INCLUDE "mpp.h"
#include "sp.h"
!-----------------------------------------------------------------------
    PARAMETER &
    (CAPA=0.28589641,RTD=57.2957795 &
    , WA=.10,WG=1.-WA,KSMUD=0)
!-------------------------CLOUD----------------------------------------
    PARAMETER &
    (A1=610.78,A2=17.2693882,A3=273.16,A4=35.86 &
!     &, PQ0=379.90516,SNOALB=0.55)
    , PQ0=379.90516)
!-------------------------CLOUD----------------------------------------
    PARAMETER &
    (IMJM=IM*JM-JM/2,JAM=6+2*(JM-10),LM1=LM-1,LP1=LM+1)

    PARAMETER &
    (SLPM=1.01325E5,EPSQ1=1.E-5,EPSQ=2.E-12,EPSO3=1.E-10,HPINC=1.E1 &
    , CLDRH0=0.80,TRESH=1.00,RNRM=1./(TRESH-CLDRH0) &
    , CLDRH2=0.90,TRESH2=1.00,RNRM2=1./(TRESH2-CLDRH2) &
    , CLAPSE=-0.0005,CLPSE=-0.0006,DCLPS=-0.0001 &
    , CM1=2937.4,CM2=4.9283,CM3=23.5518,EPS=0.622,PBOT=10000.0 &
    , STBOL=5.67E-8,PI2=2.*3.14159265,RLAG=14.8125)

    PARAMETER &
    (NB=12)
    REAL, PARAMETER :: US=1., CLIMIT=1.0E-20
!-----------------------------------------------------------------------
    PARAMETER &
    (K15=SELECTED_REAL_KIND(15))

    REAL(K15) :: PROD,DDX,EEX
!-----------------------------------------------------------------------
    LOGICAL :: &
    RUN,FIRST,RESTRT,SIGMA,CALL1,SHORT,LONG &
    ,BCLD(IDIM1:IDIM2),BTEMP1(IDIM1:IDIM2) &
    ,BITX,BITY,BITZ,BITW,BIT1,BIT2,BITC,BITS,BITCP1,BITSP1
!-------------------------CONVECTION------------------------------------
    LOGICAL :: &
    CNCLD
!-------------------------CONVECTION------------------------------------
    INCLUDE "COMM_CTLBLK.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_LOOPS.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_MASKS.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_DYNAMD.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_PHYS.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_VRBLS.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_PVRBLS.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_SOIL.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_CLDWTR.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_CNVCLD.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_INDX.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_ACMCLD.f90"
!----------------------------------------------------------------------
    INCLUDE "COMM_ACMRDL.f90"
!----------------------------------------------------------------------
    INCLUDE "COMM_ACMRDS.f90"
!-----------------------------------------------------------------------
    COMMON /CUINIT/ CURAD
    LOGICAL :: CURAD
!-----------------------------------------------------------------------
    COMMON &
    /SWRSAV/ABCFF(NB),PWTS(NB),CFCO2,CFO3,REFLO3,RRAYAV
    COMMON &
    /RD1TIM/K400,CTHK(3),LTOP(3),PTOPC(4),TAUCV(3),R1 &
    , LVL(IDIM1:IDIM2,JDIM1:JDIM2)
!-----------------------------------------------------------------------
    DIMENSION &
    TENDK (LM),CLDAMT(0:LM) &
    , PSFC  (IDIM1:IDIM2),TSKN (IDIM1:IDIM2) &
!     &, ALBEDO(IDIM1:IDIM2),XLAT(IDIM1:IDIM2),COSZ (IDIM1:IDIM2)
    , ALBDO (IDIM1:IDIM2),XLAT(IDIM1:IDIM2),COSZ (IDIM1:IDIM2) &
    , CLDCFR(IDIM1:IDIM2,3),MBOT(IDIM1:IDIM2,3) &
    , CLDF  (IDIM1:IDIM2,LP1),SLMSK (IDIM1:IDIM2) &
    , TENDS (IDIM1:IDIM2,LM),TENDL (IDIM1:IDIM2,LM) &

    , PMID  (IDIM1:IDIM2,LM),TMID  (IDIM1:IDIM2,LM) &
    , QMID  (IDIM1:IDIM2,LM),THMID(IDIM1:IDIM2,LM) &
    , OZN   (IDIM1:IDIM2,LM),POZN  (IDIM1:IDIM2,LM) &
    , MTOP(IDIM1:IDIM2,3),ICVB(IDIM1:IDIM2), ICVT(IDIM1:IDIM2) &
    , CV(IDIM1:IDIM2),SV(IDIM1:IDIM2) &

    , FLWUP (IDIM1:IDIM2),FSWDN (IDIM1:IDIM2),FSWUP (IDIM1:IDIM2) &
    , FSWDNS(IDIM1:IDIM2),FSWUPS(IDIM1:IDIM2) &
    , FLWDNS(IDIM1:IDIM2),FLWUPS(IDIM1:IDIM2) &
    , PDSL  (IDIM1:IDIM2,JDIM1:JDIM2) &
    , FNE(IDIM1:IDIM2,JDIM1:JDIM2),FSE(IDIM1:IDIM2,JDIM1:JDIM2) &
    , TL (IDIM1:IDIM2,JDIM1:JDIM2)
    DIMENSION &
    PBOTL(IDIM1:IDIM2,JDIM1:JDIM2), PTOPL(IDIM1:IDIM2,JDIM1:JDIM2) &
    ,PBOTM(IDIM1:IDIM2,JDIM1:JDIM2), PTOPM(IDIM1:IDIM2,JDIM1:JDIM2) &
    ,PBOTH(IDIM1:IDIM2,JDIM1:JDIM2), PTOPH(IDIM1:IDIM2,JDIM1:JDIM2) &
    ,TOT  (IDIM1:IDIM2,JDIM1:JDIM2)
    DIMENSION &
    CC(9),PPT(9)
    DIMENSION &
    PINT(IDIM1:IDIM2,LP1),PHALF(LP1),CSTR(IDIM1:IDIM2) &
    , EMIS(IDIM1:IDIM2,LP1), TAUC(IDIM1:IDIM2) &
    , CVB(IDIM1:IDIM2),CVT(IDIM1:IDIM2),TAUDAR(IDIM1:IDIM2)
    DIMENSION &
    CAMT(IDIM1:IDIM2,LP1),NCLDS(IDIM1:IDIM2) &
    ,ITYP(IDIM1:IDIM2,LP1),KTOP(IDIM1:IDIM2,LP1) &
    ,KBTM(IDIM1:IDIM2,LP1),RRCL(IDIM1:IDIM2,NB,LP1) &
    ,TTCL(IDIM1:IDIM2,NB,LP1),KCLD(IDIM1:IDIM2)
!--------------------CLOUD----------------------------------------------
    DIMENSION &
    CCR(IDIM1:IDIM2,LM),IW(IDIM1:IDIM2,LM),CSMID(IDIM1:IDIM2,LM) &
    ,WMID(IDIM1:IDIM2,LM),HMID(IDIM1:IDIM2,LM) &
    ,BMID(IDIM1:IDIM2),UMID(IDIM1:IDIM2) &
    ,CCMID(IDIM1:IDIM2,LM)
!--------------------CLOUD----------------------------------------------
    DATA &
    PLOMD/64200./,PMDHI/35000./,PHITP/15000./,P400/40000./ &
    , PLBTM/105000./
    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./
!----------------------------------------------------------------------
    UTIM=1.
    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=MYJS,MYJE
        DO I=MYIS,MYIE
            PDSL(I,J)=RES(I,J)*PD(I,J)
        ENDDO
    ENDDO
!**********************************************************************
!***  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) SHORT= .TRUE. 
    IF (MOD(NTSD,NRADL) == 1) LONG= .TRUE. 
    IF (MYPE == 0) THEN
        IF (SHORT) THEN
            WRITE(0,"(a)") 'RADTN: CALCULATE SHORTWAVE'
            WRITE(6,"(a)") 'RADTN: CALCULATE SHORTWAVE'
        ENDIF
        IF (LONG) THEN
            WRITE(0,"(a)") 'RADTN: CALCULATE LONGWAVE'
            WRITE(6,"(a)") 'RADTN: CALCULATE LONGWAVE'
        ENDIF
    ENDIF
    ITIMSW=0
    ITIMLW=0
    IF(SHORT)ITIMSW=1
    IF(LONG) ITIMLW=1
!-----------------------------------------------------------------------
!***
!***  FLAG FOR RESETTING CUPPT,HTOP,HBOT IN CHKOUT
!***
    IF (MOD(NTSD,NRADPP) == 1) CURAD= .TRUE. 
!***
!***  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+0.50)
    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)
    IF(SHORT)THEN
    ! omp parallel do private(i,j)
        DO J=MYJS,MYJE
            DO I=MYIS,MYIE
                CZMEAN(I,J)=0.
                TOT(I,J)=0.
            ENDDO
        ENDDO
    
        DO II=0,NRADS,NPHS
            TIMES=(NTSD-1)*DT+II*DT
            CALL ZENITH(TIMES,DAYI,HOUR)
        ! omp parallel do private(i,j)
            DO J=MYJS,MYJE
                DO I=MYIS,MYIE
                    IF(CZEN(I,J) > 0.)THEN
                        CZMEAN(I,J)=CZMEAN(I,J)+CZEN(I,J)
                        TOT(I,J)=TOT(I,J)+1.
                    ENDIF
                ENDDO
            ENDDO
        ENDDO
    ! omp parallel do private(i,j)
        DO J=MYJS,MYJE
            DO I=MYIS,MYIE
                IF(TOT(I,J) > 0.)CZMEAN(I,J)=CZMEAN(I,J)/TOT(I,J)
            ENDDO
        ENDDO
    ENDIF

! 345678901234567890123456789012345678901234567890123456789012345678901
! omp parallel do
!!$omp& private (aa,albedo,apes,bb,bcld,bit1,bit2,bitc,bitcp1,bits,
! omp& private (aa,albdo,apes,bb,bcld,bit1,bit2,bitc,bitcp1,bits,
! omp&          bitsp1,bitw,bitx,bity,bitz,bmid,btemp1,camt,
! omp&          cc1,cc2,ccmid,ccr,cfravg,cl1,cl2,
! omp&          cldamt,cldcfr,cldmax,clpfil,cosz,cr1,
! omp&          csmid,cstr,cv,cwmkl,dd,ddp,delp,denom,
! omp&          dpcl,dthdp,ee,emis,exner,fctra,
! omp&          fctrb,ff,fiq,fiw,flwdns,flwup,flwups,fswdn,fswdns,
! omp&          fswup,fswups,gg,hh,hmid,i,icvb,icvt,ir,ityp,
! omp&          iw,iwkl,j,kbt1,kbt2,kbtm,kcld,kntlyr,
! omp&          kth1,kth2,ktop,ktop1,l,l400,lbase,lin,ll,
! omp&          llbot,lltop,lml,ltrop,lvlij,malvl,mbot,mtop,
! omp&          n,nband,nbtm,nc,ncld,nclds,nktp,nlvl,
! omp&          nmod,ozn,p1,p2,pdslij,pint,
! omp&          pmid,pmod,pozn,pp,prs1,prs2,psfc,qc,qi,qint,qkl,
! omp&          qmid,qsum,qw,rqkl,rrcl,slmsk,snofac,sv,
! omp&          tauc,taudar,tcld,tendl,tends,thmid,
! omp&          tkl,tmid,tmt0,tmt15,tskn,ttcl,
! omp&          u00kl,umid,wmid,xlat)

!********************************************************************
!***  THIS IS THE BEGINNING OF THE PRIMARY LOOP THROUGH THE DOMAIN
!********************************************************************
!                        *********************
    DO 700 J = MYJS, MYJE
    !                        *********************
    
        DO 125 L=1,LM
            DO I=MYIS,MYIE
                IR=IRAD(I)
                TMID(I,L)=T(I,J,1)
                QMID(I,L)=EPSQ
                CSMID(I,L)=0.
                WMID(I,L)=0.
                CCMID(I,L)=0.
                IW(I,L)=0.
                CCR(I,L)=0.
                HMID(I,L)=0.
                OZN(I,L)=EPSO3
                TENDS(I,L)=0.
                TENDL(I,L)=0.
            ENDDO
        125 END DO
    
        DO 140 N=1,3
            DO I=MYIS,MYIE
                CLDCFR(I,N)=0.
                MTOP(I,N)=0
                MBOT(I,N)=0
            ENDDO
        140 END DO
    !***
    !***  FILL IN WORKING ARRAYS WHERE VALUES AT L=LM ARE THOSE THAT
    !***  ARE ACTUALLY AT ETA LEVEL L=LMH.
    !***
        DO 200 I=MYIS,MYIE
            IR=IRAD(I)
            LML=LMH(I,J)
            LVLIJ=LVL(I,J)
        ! lb  BMID(I)=HBM2(IR,J)
            BMID(I)=HBM2(I,J)
            UMID(I)=U00(I,J)
        
            DO L=1,LML
                PMID(I,L+LVLIJ)=AETA(L)*PDSL(I,J)+PT
                PINT(I,L+LVLIJ+1)=ETAD(L+1)*PDSL(I,J)+PT
                EXNER=(1.E5/PMID(I,L+LVLIJ))**CAPA
                TMID(I,L+LVLIJ)=T(I,J,L)
                THMID(I,L+LVLIJ)=T(I,J,L)*EXNER
                QMID(I,L+LVLIJ)=Q(I,J,L)
                WMID(I,L+LVLIJ)=CWM(I,J,L)
                HMID(I,L+LVLIJ)=HTM(I,J,L)
            ENDDO
        !***
        !***  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(LVLIJ > 0)THEN
                KNTLYR=0
            
                DO L=LVLIJ,1,-1
                    KNTLYR=KNTLYR+1
                    PMID(I,L)=PT-REAL(2*KNTLYR-1)*0.5*HPINC
                    PINT(I,L+1)=PMID(I,L)+0.5*HPINC
                    EXNER=(1.E5/PMID(I,L))**CAPA
                    THMID(I,L)=TMID(I,L)*EXNER
                ENDDO
            ENDIF
        
            IF(LVLIJ == 0) THEN
                PINT(I,1)=PT
            ELSE
                PINT(I,1)=PMID(I,1)-0.5*HPINC
            ENDIF
        200 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 250 I=MYIS,MYIE
            PSFC(I)=PD(I,J)+PT
            APES=(PSFC(I)*1.E-5)**CAPA
            TSKN(I)=THS(I,J)*APES*(1.-2.*SM(I,J))
            SLMSK(I)=SM(I,J)
        
        ! ----------------------------------------------------------------------
        ! turn off snow albedo calculation since it is now calculated in SFLX.
        !      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))
            ALBDO(I)=ALBEDO(I,J)
        
            XLAT(I)=GLAT(I,J)*RTD
            COSZ(I)=CZMEAN(I,J)
        250 END DO
    !-----------------------------------------------------------------------
    !*******************STRATIFORM CLOUD SECTION***************************
    !-----------------------------------------------------------------------
    !  CALCULATE STRATIFORM CLOUD COVERAGE AT EACH MODEL GRID POINT WHICH
    !  WILL BE USED IN THE MODEL RADIATION PARAMETERIZATION SCHEME.
    !-----------------------------------------------------------------------
    !------------------QW, QI AND QINT--------------------------------------
        DO 280 I=MYIS,MYIE
            LML=LMH(I,J)
            LVLIJ=LVL(I,J)
        
            DO 275 L=1,LML
                LL=L+LVLIJ
                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
                ENDIF
            
                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 > 0. .OR. CWMKL > CLIMIT)THEN
                        IW(I,LL)=1
                    ELSE
                        IW(I,LL)=0
                    ENDIF
                ENDIF
            
                IF(TMT0 >= 0.)THEN
                    IW(I,LL)=0
                ENDIF
            
                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
                ENDIF
            
                IWKL=IW(I,LL)
            
            !----------------THE SATURATION SPECIFIC HUMIDITY------------------------
            
                FIW=FLOAT(IWKL)
                QC=(1.-FIW)*QINT+FIW*QI
            
            !----------------THE RELATIVE HUMIDITY----------------------------------
            
                IF(QC <= EPSQ1 .OR. QKL <= EPSQ1)THEN
                    RQKL=0.
                ELSE
                    RQKL=QKL/QC
                ENDIF
            
            !----------------CLOUD COVER RATIO CCR----------------------------------
            
                IF(RQKL >= 0.9999)THEN
                    CCR(I,LL)=AMIN1(US,RQKL)
                
                !--- Ferrier, 6/17/02:  Emergency change to eliminate very small cloud
                !    amounts (soon to be replaced by corrections from ETAZ parallel)
                
                ELSE IF (CWMKL > 1.E-7) THEN
                    ARG=-1000.*CWMKL/(US-RQKL)
                    ARG=AMAX1(ARG,-25.)
                    CCR_dum=RQKL*(1.-EXP(ARG))
                    IF (CCR_dum > .01) CCR(I,LL)=CCR_dum
                ENDIF
                CSMID(I,LL)=AMIN1(US,CCR(I,LL))
            !----------------------------------------------------------------------
            275 END DO
        280 END DO
    !----------------------------------------------------------------------
    !**********************************************************************
    ! NOW CHECK THE CLOUDS PRODUCED ABOVE TO MAKE SURE THEY ARE GOOD
    ! ENOUGH FOR RADIATION CALCULATIONS
    !**********************************************************************
    !***
    !***  NO STRATIFORM CLOUDS FOR THIS TYPE
    !***
        DO 350 I=MYIS,MYIE
        
            LML=LMH(I,J)
            LVLIJ=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+LVLIJ
                    DDP=PSFC(I)-PMID(I,LL)
                    IF(DDP >= PBOT) GO TO 290
                    CSMID(I,LL)=0.
                ENDDO
                290 CONTINUE
            ENDIF
        !***
        !***  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+LVLIJ
                    DTHDP=(THMID(I,LL-1)-THMID(I,LL)) &
                    /(PMID(I,LL-1)-PMID(I,LL))
                    IF(DTHDP <= CLAPSE)THEN
                        LBASE=LL
                        GO TO 300
                    ENDIF
                ENDDO
                300 CONTINUE
            
            !***  CHECK 2 LAYERS ABOVE LBASE FOR DRYNESS
            
                IF(CSMID(I,LBASE-1) <= 0. .AND. CSMID(I,LBASE-2) <= 0. &
                 .AND. LBASE < LM)THEN
                    IF(DTHDP > CLPSE)THEN
                        CLPFIL=1.-((CLPSE-DTHDP)/DCLPS)
                    ELSE
                        CLPFIL=1.
                    ENDIF
                
                    DO L=1,LML
                        LL=LML-L+1+LVLIJ
                        DDP=PSFC(I)-PMID(I,LL)
                        IF(DDP >= PBOT) GO TO 310
                        CSMID(I,LL)=CSMID(I,LL)*CLPFIL
                    ENDDO
                    310 CONTINUE
                
                !***  IF NO INVERSION OR IF CLDS EXIST IN EITHER OF THE 2 LAYERS ABOVE
                !***  INVERSION, ZERO OUT CLOUD BELOW PBOT
                
                ELSE
                    DO L=1,LML
                        LL=LML-L+1+LVLIJ
                        DDP=PSFC(I)-PMID(I,LL)
                        IF(DDP >= PBOT) GO TO 320
                        CSMID(I,LL)=0.
                    ENDDO
                    320 CONTINUE
                
                ENDIF
            !------------
            ENDIF
        !------------
        !***
        !***  REMOVE HIGH CLOUDS ABOVE THE TROPOPAUSE
        !***
            L400=LM
            DO L=1,LML
                LL=LML-L+1+LVLIJ
                IF(PMID(I,LL) <= 40000.0)THEN
                    L400=LL
                    GO TO 330
                ENDIF
            ENDDO
            330 CONTINUE
        
            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
                    GOTO 340
                ENDIF
            ENDDO
            340 IF(LTROP < LM)THEN
                DO LL=LTROP,1,-1
                    CSMID(I,LL)=0.
                ENDDO
            ENDIF
        350 END DO
    
    !*********************************************************************
    !*****************END OF STRATIFORM CLOUD SECTION*****************
    !------------------------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
    
    !***
    !*** NON-PRECIPITATING CLOUD FRACTION OF 20 PERCENT IS ADDED AT
    !*** AT POINTS WHERE THE SHALLOW AND DEEP CONVECTIONS ACCUR.
    !***                                        Q. ZHAO  97-5-2
    
    !   COMPUTE THE CONVECTIVE CLOUD COVER FOR RADIATION
    
    !-----------------------------------------------------------------
        IF(CNCLD)THEN
        !-----------------------------------------------------------------
            DO 375 I=MYIS,MYIE
                IF(HBOT(I,J)-HTOP(I,J) > 1.0)THEN
                !       IF(HTOP(I,J).LT.HBOT(I,J))THEN
                    SV(I)=0.0
                ELSE
                    SV(I)=0.0
                ENDIF
            
                PMOD=CUPPT(I,J)*24.0*1000.0/CLSTP
                NMOD=0
            
                DO NC=1,9
                    IF(PMOD > PPT(NC)) NMOD=NC
                ENDDO
            
            !***  CLOUD TOPS AND BOTTOMS COME FROM CUCNVC
            !***  ADD LVL TO BE CONSISTENT WITH OTHER WORKING ARRAYS
            
                IF(NMOD == 0)THEN
                    CV(I)=0.
                ELSEIF(NMOD == 9)THEN
                    CV(I)=CC(9)
                ELSE
                    CC1=CC(NMOD)
                    CC2=CC(NMOD+1)
                    P1=PPT(NMOD)
                    P2=PPT(NMOD+1)
                    CV(I)=CC1+(CC2-CC1)*(PMOD-P1)/(P2-P1)
                ENDIF
            
                CV(I)=AMAX1(SV(I),CV(I))
                CV(I)=AMIN1(1.0,CV(I))
            
                IF(CV(I) == 0.0)THEN
                    ICVT(I)=0
                    ICVB(I)=0
                ELSE
                    ICVT(I)=INT(HTOP(I,J)+0.50)+LVL(I,J)
                    ICVB(I)=INT(HBOT(I,J)+0.50)+LVL(I,J)
                ENDIF
            375 END DO
        !***
        !***  MAKE SURE CLOUDS ARE DEEP ENOUGH
        !***
            DO I=MYIS,MYIE
                BCLD(I)=CV(I) > 0. .AND. &
                (ICVB(I)-ICVT(I)) >= 1
                BTEMP1(I)=BCLD(I)
            ENDDO
        !***
        !*** COMPUTE CONVECTIVE CLOUD FRACTION
        !***
            DO 390 I=MYIS,MYIE
                IF(BCLD(I)) THEN
                    LML=LMH(I,J)
                    LVLIJ=LVL(I,J)
                
                    DO L=1,LML
                        LL=L+LVLIJ
                        IF(LL > ICVB(I) .OR. LL < ICVT(I))THEN
                            CCMID(I,LL)=0.
                        ELSE
                            CCMID(I,LL)=CV(I)
                        ENDIF
                        CCMID(I,LL)=AMIN1(1.0,CCMID(I,LL))
                    ENDDO
                ENDIF
            390 END DO
        !***
        !***  REMOVE HIGH CLOUDS ABOVE THE TROPOPAUSE
        !***
            L400=LM
            DO 425 I=MYIS,MYIE
                LML=LMH(I,J)
                LVLIJ=LVL(I,J)
            
                DO L = 1, LML
                    LL=LML-L+1+LVLIJ
                    IF(PMID(I,LL) <= 40000.0)THEN
                        L400=LL
                        GO TO 400
                    ENDIF
                ENDDO
                400 CONTINUE
            
                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
                        GOTO 410
                    ENDIF
                ENDDO
                410 IF(LTROP < LM)THEN
                    DO LL=LTROP,1,-1
                        CCMID(I,LL)=0.
                    ENDDO
                ENDIF
            425 END DO
        !***
        !-----------------------------------------------------------------
        ENDIF
    !------------------------CONVECTION--------------------------------
    !*****************END OF CONVECTIVE CLOUD SECTION*****************
    !*********************************************************************
    !***
    !***  DETERMINE THE FRACTIONAL CLOUD COVERAGE FOR HIGH, MID
    !***  AND LOW OF CLOUDS FROM THE CLOUD COVERAGE AT EACH LEVEL
    !***
    !***  NOTE: THIS IS FOR DIAGNOSTICS ONLY!!!
    !***
    !***
        DO 500 I=MYIS,MYIE
        
            CSTR(I)=0.0
        
            DO L=0,LM
                CLDAMT(L)=0.
            ENDDO
        
        !***  NOW GOES LOW, MIDDLE, HIGH
        
            DO 480 NLVL=1,3
                CLDMAX=0.
                MALVL=LM
                LLTOP=LTOP(NLVL)+LVL(I,J)
            !***
            !***  GO TO THE NEXT CLOUD LAYER IF THE TOP OF THE CLOUD-TYPE IN
            !***  QUESTION IS BELOW GROUND OR IS IN THE LOWEST LAYER ABOVE GROUND.
            !***
                IF(LLTOP >= LM)GO TO 480
            
                IF(NLVL > 1)THEN
                    LLBOT=LTOP(NLVL-1)-1+LVL(I,J)
                    LLBOT=MIN(LLBOT,LM1)
                ELSE
                    LLBOT=LM1
                ENDIF
            
                DO 435 L=LLTOP,LLBOT
                    CLDAMT(L)=AMAX1(CSMID(I,L),CCMID(I,L))
                    IF(CLDAMT(L) > CLDMAX)THEN
                        MALVL=L
                        CLDMAX=CLDAMT(L)
                    ENDIF
                435 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 450 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)GO TO 450
                !***
                !***  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
                    ENDIF
                
                    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
                        ENDIF
                    ENDIF
                    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)GO TO 450
                
                    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
                    ENDIF
                
                    IF(BITW)THEN
                        KBT1=KBT2
                        KTH1=KTH2
                        CL1=CL2
                    ENDIF
                
                    IF(BITY)THEN
                        KBT2=LLBOT
                        KTH2=0
                        CL2=0.0
                    ENDIF
                450 END DO
            !***
                CLDCFR(I,NLVL)=AMIN1(1.0,CL1)
                MTOP(I,NLVL)=MIN(KBT1,KBT1-KTH1+1)
                MBOT(I,NLVL)=KBT1
            480 END DO
        500 END DO
    !***
    !***  SET THE UN-NEEDED TAUDAR TO ONE
    !***
        DO I=MYIS,MYIE
            TAUDAR(I)=1.0
        ENDDO
    !----------------------------------------------------------------------
    ! 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 600 I=MYIS,MYIE
            LML=LMH(I,J)
            LVLIJ=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
            ENDDO
        
            DO 510 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
                ENDDO
            510 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 540 L=2,LML
                LL=LML-L+1+LVLIJ
                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
                        ENDIF
                    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
                            ENDIF
                        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))
                            ENDIF
                        ENDIF
                    ENDIF
                !-------------------
                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
                    ENDIF
                !-------------------
                ENDIF
            !-------------------
            540 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 580 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 560 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)
                                ENDIF
                            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
                                ENDIF
                            ENDIF
                        ENDIF
                    560 END DO
                
                    IF(BITX)EMIS(I,NC)=1.0-EXP(-0.75*TAUC(I))
                    IF(QSUM >= EPSQ1)THEN
                    
                        DO 570 NBAND=1,NB
                            IF(BITX)THEN
                                PROD=ABCFF(NBAND)*QSUM
                                DDX=TAUC(I)/(TAUC(I)+PROD)
                                EEX=1.0-DDX
                                IF(ABS(EEX) >= 1.E-8)THEN
                                    DD=DDX
                                    EE=EEX
                                    FF=1.0-DD*0.85
                                    AA=MIN(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)=MAX(0.1E-5,(BB-1.0)*(1.0-AA*AA)/DD)
                                    TTCL(I,NBAND,NC)=AMAX1(0.1E-5,4.0*GG*AA/DD)
                                ENDIF
                            ENDIF
                        570 END DO
                    ENDIF
                580 END DO
            
            ENDIF
        
        600 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=MYIS,MYIE
                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
            ENDDO
        ENDDO
    
        CALL OZON2D(LM,POZN,XLAT,RSIN1,RCOS1,RCOS2,OZN)

    !***
    !***  NOW THE VARIABLES REQUIRED BY RADFS HAVE BEEN CALCULATED.
    !***
    !----------------------------------------------------------------------
    !***
    !***  CALL THE GFDL RADIATION DRIVER
    !***
    !***
        CALL RADFS &
    !     1     (PSFC,PMID,PINT,QMID,TMID,OZN,TSKN,SLMSK,ALBEDO,XLAT
        (PSFC,PMID,PINT,QMID,TMID,OZN,TSKN,SLMSK,ALBDO,XLAT &
        ,     CAMT,ITYP,KTOP,KBTM,NCLDS,EMIS,RRCL,TTCL &
        ,     COSZ,TAUDAR,1 &
        ,     1,0 &
        ,     ETAD,AETA,ITIMSW,ITIMLW,JD,R1,HOUR,TENDS,TENDL &
        ,     FLWUP,FSWUP,FSWDN,FSWDNS,FSWUPS,FLWDNS,FLWUPS)
    !----------------------------------------------------------------------
        DO 650 I=MYIS,MYIE
            PDSLIJ=PDSL(I,J)
            PMOD=CUPPT(I,J)*24.0*1000.0/CLSTP
            CFRACL(I,J)=CLDCFR(I,1)
            CFRACM(I,J)=CLDCFR(I,2)
            CFRACH(I,J)=CLDCFR(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)))
        
        !--- Prevent double-counting statistics when restarting model
        
            IF (NTSD == 1 .OR. NTSD > NSTART+1) THEN
                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
                    ENDIF
                ELSE
                    ACFRST(I,J)=ACFRST(I,J)+CFRAVG
                    NCFRST(I,J)=NCFRST(I,J)+1
                ENDIF
            ENDIF
        650 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 660 I=MYIS,MYIE
            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)GO TO 660
            ENDDO
        660 END DO
    !***
    !***  SUM THE LW INCOMING AND SW RADIATION (W/M**2) FOR RADIN.
    !***
        DO 675 I=MYIS,MYIE
            IF(LONG)THEN
                SIGT4(I,J)=STBOL*TMID(I,LM)*TMID(I,LM)* &
                TMID(I,LM)*TMID(I,LM)
            ENDIF
        
        !***  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)
!if (mype == 2) write(*,*) "TESTERADT", FLWDNS(I),FLWUPS(I),FLWUP(I)
            ENDIF
            IF(SHORT)THEN
                RSWIN(I,J) =FSWDNS(I)
                RSWOUT(I,J)=FSWUPS(I)
                RSWTOA(I,J)=FSWUP(I)
!if (mype == 2) write(*,*) "TESTERADN", FSWDNS(I),FSWUPS(I),FSWUP(I)
            ENDIF
        675 END DO
    !***
    !***  THIS ROW IS FINISHED. GO TO NEXT
    !***
    !                        *********************
    700 END DO
!                        *********************
!----------------------------------------------------------------------
!***
!***  CALLS TO RADIATION THIS TIME STEP ARE COMPLETE.
!***
!----------------------------------------------------------------------
!----------------------------------------------------------------------
!***
!***  HORIZONTAL SMOOTHING OF TEMPERATURE TENDENCIES
!***
!----------------------------------------------------------------------
    IF(SHORT) THEN
        DO 800 L=1,LM
            CALL ZERO2(TL)
            CALL ZERO2(FNE)
            CALL ZERO2(FSE)
        
            IF(KSMUD >= 1)THEN
                DO 750 KS=1,KSMUD
                
                    DO J=MYJS,MYJE
                        DO I=MYIS,MYIE
                            TL(I,J)=RSWTT(I,J,L)*HTM(I,J,L)
                        ENDDO
                    ENDDO
                
                    DO J=MYJS,MYJE
                        DO I=MYIS,MYIE
                            FNE(I,J)=(TL(I+IHE(J),J+1)-TL(I,J)) &
                            *HTM(I,J,L)*HTM(I+IHE(J),J+1,L)
                        ENDDO
                    ENDDO
                
                    DO J=MYJS1,MYJE
                        DO I=MYIS,MYIE
                            FSE(I,J)=(TL(I+IHE(J),J-1)-TL(I,J)) &
                            *HTM(I+IHE(J),J-1,L)*HTM(I,J,L)
                        ENDDO
                    ENDDO
                
                    DO J=MYJS2,MYJE2
                        DO I=MYIS,MYIE
                            TL(I,J)=(FNE(I,J)-FNE(I+IHW(J),J-1) &
                            +FSE(I,J)-FSE(I+IHW(J),J+1)) &
                            *HBM2(I,J)*0.125+TL(I,J)
                        ENDDO
                    ENDDO
                
                    DO J=MYJS,MYJE
                        DO I=MYIS,MYIE
                            RSWTT(I,J,L)=TL(I,J)
                        ENDDO
                    ENDDO
                
                750 END DO
            ENDIF
        
        800 END DO
    ENDIF
!----------------------------------------------------------------------

    IF(LONG)THEN
    
        DO 900 L=1,LM
            CALL ZERO2(TL)
            CALL ZERO2(FNE)
            CALL ZERO2(FSE)
        
            IF(KSMUD >= 1)THEN
                DO 850 KS=1,KSMUD
                
                    DO J=MYJS,MYJE
                        DO I=MYIS,MYIE
                            TL(I,J)=RLWTT(I,J,L)*HTM(I,J,L)
                        ENDDO
                    ENDDO
                
                    DO J=MYJS,MYJE1
                        DO I=MYIS,MYIE
                            FNE(I,J)=(TL(I+IHE(J),J+1)-TL(I,J)) &
                            *HTM(I,J,L)*HTM(I+IHE(J),J+1,L)
                        ENDDO
                    ENDDO
                
                    DO J=MYJS1,MYJE
                        DO I=MYIS,MYIE
                            FSE(I,J)=(TL(I+IHE(J),J-1)-TL(I,J)) &
                            *HTM(I+IHE(J),J-1,L)*HTM(I,J,L)
                        ENDDO
                    ENDDO
                
                    DO J=MYJS2,MYJE2
                        DO I=MYIS,MYIE
                            TL(I,J)=(FNE(I,J)-FNE(I+IHW(J),J-1) &
                            +FSE(I,J)-FSE(I+IHW(J),J+1)) &
                            *HBM2(I,J)*0.125+TL(I,J)
                        ENDDO
                    ENDDO
                
                    DO J=MYJS,MYJE
                        DO I=MYIS,MYIE
                            RLWTT(I,J,L)=TL(I,J)
                        ENDDO
                    ENDDO
                
                850 END DO
            ENDIF
        900 END DO
    ENDIF
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
    RETURN
    END SUBROUTINE RADTN
