radlsw.F90

Go to the documentation of this file.

SUBROUTINE radlsw &
      &( kidia, kfdia , klon , ktdia, klev  , kmode, kaer, kbox, nbox &
      &, ndump, klwrad &
      &, prii0 &
      &, paer , palbd , palbp, paph , pap &
      &, pcco2, pfrcl , pdp  , pemis, pemiw , plsm , pmu0, pozon &
      &, pq   , pqiwp , pqlwp, psqiw, psqlw , pqs  , pqrain, praint &
      &, prlvri,prlvrl, pth  , pt   , pts   , pnbas, pntop &
      &, pemit, pfct  , pflt , pfcs , pfls  , pfrsod, psudu, puvdf, pparf &
      &, pfdct, pfuct , pfdlt, pfult, pfdcs , pfucs , pfdls, pfuls &
      &, ztau , ztauint &
      &, aswbox, olrbox, slwbox, sswbox, taubox, pclbx &
! #DB &, k2iii, k2jjj &
      &)

!**** *RADLSW* - INTERFACE TO ECMWF LW AND SW RADIATION SCHEMES

!     PURPOSE.
!     --------
!           CONTROLS RADIATION COMPUTATIONS

!**   INTERFACE.
!     ----------

!        EXPLICIT ARGUMENTS :
!        --------------------
! PAER   : (KLON,6,KLEV)     ; OPTICAL THICKNESS OF THE AEROSOLS
! PALBD  : (KLON,NSW)        ; SURF. SW ALBEDO FOR DIFFUSE RADIATION
! PALBP  : (KLON,NSW)        ; SURF. SW ALBEDO FOR PARALLEL RADIATION
! PAPH   : (KLON,KLEV+1)     ; HALF LEVEL PRESSURE
! PAP    : (KLON,KLEV)       ; FULL LEVEL PRESSURE
! PCCO2  :                   ; CONCENTRATION IN CO2 (PA/PA)
! PFRCL  : (KLON,KLEV)       ; CLOUD FRACTIONAL COVER
! PDP    : (KLON,KLEV)       ; LAYER PRESSURE THICKNESS
! PEMIS  : (KLON)            ; SURFACE LW EMISSIVITY
! PEMIW  : (KLON)            ; SURFACE LW WINDOW EMISSIVITY
! PLSM   : (KLON)            ; LAND-SEA MASK
! PMU0   : (KLON)            ; SOLAR ANGLE
! PNBAS  : (KLON)            ; INDEX OF BASE OF CONVECTIVE LAYER
! PNTOP  : (KLON)            ; INDEX OF TOP OF CONVECTIVE LAYER
! POZON  : (KLON,KLEV)       ; CONCENTRATION IN OZONE (PA/PA)
! PQ     : (KLON,KLEV)       ; SPECIFIC HUMIDITY PA/PA
! PQIWP  : (KLON,KLEV)       ; SOLID  WATER KG/KG
! PQLWP  : (KLON,KLEV)       ; LIQUID WATER KG/KG
! PQS    : (KLON,KLEV)       ; SATURATION WATER VAPOR  KG/KG
! PQRAIN : (KLON,KLEV)       ; RAIN WATER KG/KG
! PRAINT : (KLON,KLEV)       ; RAIN RATE (m/s)
! PRLVRI : (KLON,KLEV)       ; RELATIVE VARIANCE OF ICE WATER
! PRLVRL : (KLON,KLEV)       ; RELATIVE VARIANCE OF LIQUID WATER
! PTH    : (KLON,KLEV+1)     ; HALF LEVEL TEMPERATURE
! PT     : (KLON,KLEV)       ; FULL LEVEL TEMPERATURE
! PTS    : (KLON)            ; SURFACE TEMPERATURE
!     ==== OUTPUTS ===
! PFCT   : (KLON,KLEV+1)     ; CLEAR-SKY LW NET FLUXES
! PFLT   : (KLON,KLEV+1)     ; TOTAL LW NET FLUXES
! PFCS   : (KLON,KLEV+1)     ; CLEAR-SKY SW NET FLUXES
! PFLS   : (KLON,KLEV+1)     ; TOTAL SW NET FLUXES
! PFRSOD : (KLON)            ; TOTAL-SKY SURFACE SW DOWNWARD FLUX
! PEMIT  : (KLON)            ; SURFACE TOTAL LONGWAVE EMISSIVITY
! PSUDU  : (KLON)            ; SOLAR RADIANCE IN SUN'S DIRECTION
! PUVDF  : (KLON)            ; SURFACE DOWNWARD U.V. RADIATION
! PPARF  : (KLON)            ; PHOTOSYNTHETICALLY ACTIVE RADIATION

!        IMPLICIT ARGUMENTS :   NONE
!        --------------------

!     METHOD.
!     -------
!        SEE DOCUMENTATION

!     EXTERNALS.
!     ----------

!     REFERENCE.
!     ----------
!        ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE IFS

!     AUTHORS.
!     --------
!        J.-J. MORCRETTE         *ECMWF*

!     MODIFICATIONS.
!     --------------
!        ORIGINAL : 88-02-04
!        J.-J. MORCRETTE 94-11-15 DIRECT/DIFFUSE SURFACE ALBEDO
!        08/96: J.-J. Morcrette/Ph. Dandin: tests of eff. radius param.
!        9909 : JJMorcrette effect.radius + inhomogeneity factors
!        JJMorcrette 990128 : sunshine duration
!        JJMorcrette : 990831 RRTM-140gp
!-----------------------------------------------------------------------

#include "tsmbkind.h"

!USE YOMCT3   , ONLY : NSTEP
USE yomcst   , ONLY : rg       ,rd       ,rtt      ,rpi
USE yoerad   , ONLY : nsw      ,lrrtm    ,linhom, &
            &loifuec, ltempds, lowasyf, lowhsss, nradip, nradlp, &
            &niceopt, nliqopt, novlp  , nhowinh, rminice
USE yoelw    , ONLY : nsil     ,ntra     ,nua      ,tstand   ,xp
USE yoesw    , ONLY : ryfwca   ,ryfwcb   ,ryfwcc   ,ryfwcd   ,&
            &ryfwce   ,ryfwcf   ,rebcua   ,rebcub   ,rebcuc   ,&
            &rebcud   ,rebcue   ,rebcuf   ,rebcui   ,rebcuj   ,&
            &rebcug   ,rebcuh   ,rhsavi   ,rfulio   ,rflaa0   ,&
            &rflaa1   ,rflbb0   ,rflbb1   ,rflbb2   ,rflbb3   ,&
            &rflcc0   ,rflcc1   ,rflcc2   ,rflcc3   ,rfldd0   ,&
            &rfldd1   ,rfldd2   ,rfldd3   ,rfuaa0   ,rfuaa1   ,&
            &rfubb0   ,rfubb1   ,rfubb2   ,rfubb3   ,rfucc0   ,&
            &rfucc1   ,rfucc2   ,rfucc3   ,rfueta   ,raswca   ,&
            &raswcb   ,raswcc   ,raswcd   ,raswce   ,raswcf   ,&
            &rlinli
USE yoerdu   , ONLY : nuaer    ,ntraer   ,replog   ,repsc    ,diff
USE yoerdi   , ONLY : repclc
USE yoethf   , ONLY : rtice
USE yoephli  , ONLY : lphylin
USE yoerrtwn , ONLY : ng        ,nspa      ,nspb      ,wavenum1  ,&
           &wavenum2  ,delwave   ,totplnk   ,totplk16
USE yoedbug  , ONLY : ldebug


IMPLICIT NONE


!     DUMMY INTEGER SCALARS
integer_m :: kaer
integer_m :: kfdia
integer_m :: kidia
integer_m :: klev
integer_m :: klon
integer_m :: kmode
integer_m :: ktdia
integer_m :: kbox
integer_m :: nbox
integer_m :: ndump, klwrad

!     DUMMY REAL SCALARS
real_b :: prii0



!     -----------------------------------------------------------------

!*       0.1   ARGUMENTS.
!              ----------
real_b :: palbd(klon,nsw) , palbp(klon,nsw)
real_b :: pemis(klon)     , pemiw(klon)
real_b :: plsm(klon)      , pmu0(klon)
real_b :: pcco2           , pozon(klon,klev)
real_b :: pts(klon)       , pnbas(klon)     , pntop(klon)
real_b :: pt(klon,klev)  , pap(klon,klev)
real_b :: pth(klon,klev+1), paph(klon,klev+1)
real_b :: pdp(klon,klev)
real_b :: pq(klon,klev)  , pqs(klon,klev)
real_b :: pqiwp(klon,klev), pqlwp(klon,klev), pqrain(klon,klev)
real_b :: praint(klon,klev)
real_b :: prlvri(klon,klev),prlvrl(klon,klev)
real_b :: psqiw(klon,klev), psqlw(klon,klev)
real_b :: pfrcl(klon,klev), pclfr(klon,klev), pclbx(klon,100,klev)
real_b :: paer(klon,6,klev)

! #DB integer :: k2iii(KLON),k2jjj(KLON),kio,kjo

!     ==== COMPUTED IN RADLSW ===
real_b :: pfcs(klon,klev+1), pfct(klon,klev+1)
real_b :: pfls(klon,klev+1), pflt(klon,klev+1)
real_b :: pfrsod(klon)     , pemit(klon)
real_b :: psudu(klon)      , puvdf(klon)        , pparf(klon)
real_b :: pfdct(klon,klev+1), pfuct(klon,klev+1)
real_b :: pfdlt(klon,klev+1), pfult(klon,klev+1)
real_b :: pfdcs(klon,klev+1), pfucs(klon,klev+1)
real_b :: pfdls(klon,klev+1), pfuls(klon,klev+1)

real_b :: aswbox(klon, 100), olrbox(klon, 100)
real_b :: slwbox(klon, 100), sswbox(klon, 100), taubox(klon, 100)

!     -----------------------------------------------------------------

!*       0.2   LOCAL ARRAYS.
!              -------------
!     -----------------------------------------------------------------

!-- ARRAYS FOR LOCAL VARIABLES -----------------------------------------

integer_m :: ibas(klon)     , itop(klon)

! #DB integer :: jkjllw, jkjlsw, JAERmin, JAERmax, jTAUCLDmin, jTAUCLDmax
! #DB real :: PAERmin, PAERmax, TAUCLDmin, TAUCLDmax

real_b ::&
    &zalbd(klon,nsw)    , zalbp(klon,nsw)&
  &, zcg(klon,nsw,klev) , zomega(klon,nsw,klev)&
  &, ztau(klon,nsw,klev) &
  &, ztaucld(klon,klev,16), ztclear(klon)
real_b ::&
    &zcldld(klon,klev)  , zcldlu(klon,klev)&
  &, zcldsw(klon,klev)  , zcld0(klon,klev)&
  &, zdt0(klon)        &
  &, zemis(klon)        , zemiw(klon)&
  &, zflux(klon,2,klev+1)                 , zfluc(klon,2,klev+1)&
  &, zfiwp(klon)        , zflwp(klon)      , zfrwp(klon)&
  &, ziwc(klon)         , zlwc(klon)&
  &, zbicfu(klon)       , zkicfu1(klon)    , zkicfu2(klon)&
!cc            , ZRWC(KLON)
  &, zmu0(klon)         , zoz(klon,klev)   , zozn(klon,klev)&
  &, zozon(klon,klev)   , zpmb(klon,klev+1), zpsol(klon)&
  &, ztave(klon,klev)  , ztl(klon,klev+1)&
  &, zview(klon)
real_b ::&
    &zfcdwn(klon,klev+1), zfcup(klon,klev+1)&
  &, zfsdwn(klon,klev+1), zfsup(klon,klev+1)&
  &, zfsupn(klon)       , zfsupv(klon)&
  &, zfcupn(klon)       , zfcupv(klon)&
  &, zfsdnn(klon)       , zfsdnv(klon)&
  &, zfcdnn(klon)       , zfcdnv(klon)&
  &, zcoolr(klon,klev)  , zcoolc(klon,klev)&
  &, zheatr(klon,klev)  , zheatc(klon,klev)
real_b ::&
    &zalfice(klon)      , zgamice(klon)     , zbice(klon),  zdesr(klon) &
  &, zradip(klon)       , zradlp(klon)      , zcfudg(klon)&
!cc           , ZRADRD(KLON)
  &, zraint(klon)       , zres(klon)&
  &, ztice(klon)        , zemit(klon)       , ztauint(klon)
real_b :: zsudu(klon)   , zuvdf(klon)       , zparf(klon),  zcol(klon) &
  &, ztcc(klon)         , ztca(klon)

!-- box-type arrays

real_b :: cpfcs(klon,klev+1) , cpfct(klon,klev+1)
real_b :: cpfls(klon,klev+1) , cpflt(klon,klev+1)
real_b :: cpfrsod(klon)      , cpemit(klon)
real_b :: cpsudu(klon)       , cpuvdf(klon)       , cpparf(klon)
real_b :: cpfdct(klon,klev+1), cpfuct(klon,klev+1)
real_b :: cpfdlt(klon,klev+1), cpfult(klon,klev+1)
real_b :: cpfdcs(klon,klev+1), cpfucs(klon,klev+1)
real_b :: cpfdls(klon,klev+1), cpfuls(klon,klev+1)

!     LOCAL INTEGER SCALARS
integer_m :: ikl, jae, jk, jkl, jklp1, jkp1, jl, jnu, jrtm, jsw &
  &, nboxl, icbox, imov, indlay

!     LOCAL LOGICAL SCALARS
LOGICAL :: LLINTRP

!     LOCAL REAL SCALARS
real_b :: zasymx, zdiffd, zgi, zgl, zgr, ziwgkg, zlwgkg,&
          &zmsaid, zmsaiu, zmsald, zmsalu, zmtconv, &
          &zmtfudg, zlwfudg, zswfudg, zmultl, zoi, zol, zomgmx, zor, &
          &zrmuz, zrwgkg, ztaud, ztaumx, ztempc, &
          &ztoi, ztol, ztor, zzfiwp, zzflwp, zdpog, zpodt
real_b :: zalnd, zasea, zd, zden, zntot, znum, zratio, zcoeff, z1radi,&
          &z1radl, zbetai, zomgi, zomgp, zfdel, zwght, zvi, zvl, zvr
real_b :: zasw, zolr, zslw, zssw, zmulti, zaiwc, zbiwc,&
          &zdice, zfsr, zlgiwc, ztcels, ztblay, zaddplk, zplanck
real_b :: ztol1, ztoi1, ztor1


!     -----------------------------------------------------------------

! #DB kio = 66
! #DB kjo = 53

!if (NDUMP.LE.3) then
!  JL=KIDIA
!  DO jk=1,klev
!    print 9104,jk,PAPH(JL,JK),PTH(JL,JK),PAP(JL,JK),PT(JL,JK)&
!    &            ,PDP(JL,JK)& 
!    &            ,PQ(JL,JK),PFRCL(JL,JK),PQIWP(JL,JK),PQLWP(JL,JK)&
!    &            ,POZON(JL,JK),PQS(JL,JK)
9104 format(1x,i3,f9.1,f8.2,f9.1,f8.2,f9.1,e10.3,f7.4,4e10.3)
!  ENDDO 
!  jk=klev+1
!  print 9104,jk,PAPH(JL,JK),PTH(JL,JK)
!  print 9105,PTS(JL),(PALBD(JL,JSW),PALBP(JL,JSW),JSW=1,NSW)
9105 FORMAT(13x,f8.2,12f8.4)
!end if

!print *,'NICEOPT, NLIQOPT, NRADIP, NRADLP',NICEOPT,NLIQOPT,NRADIP,NRADLP

!-- compute total cloud cover
DO jl=kidia,kfdia
  ztcc(jl)=1.-pfrcl(jl,1)
  ztca(jl)=0.
END DO
DO jk=2,klev
  DO jl=kidia,kfdia  
    ztcc(jl)=ztcc(jl)*(1.-max(pfrcl(jl,jk),pfrcl(jl,jk-1))) &
    & /(1.-min(pfrcl(jl,jk-1),1.-repclc))
  END DO
END DO
DO jl=kidia,kfdia
  ztcc(jl)=1.-ztcc(jl)
END DO

!JL=KIDIA
!print 9106,ZTCC(JL)
9106 format(1x,'TCC :',f7.4)
!print 9107,LINHOM,NHOWINH
9107 format(1x,'LINHOM=',l8,' NHOWINH=',i2)
  





!*         1.     SET-UP INPUT QUANTITIES FOR RADIATION
!                 -------------------------------------

IF (.NOT.linhom) THEN
  zmtfudg=1.0_jprb
  zmtconv=1.0_jprb
  zswfudg=1.0_jprb
  zlwfudg=1.0_jprb
ELSE IF (linhom) THEN
  IF (nhowinh.EQ.1) THEN  
    zmtfudg=0.7_jprb
    zmtconv=0.7_jprb
    zswfudg=0.7_jprb
    zlwfudg=0.7_jprb    
  ELSE
    zmtfudg=1.0_jprb
    zmtconv=1.0_jprb
    zswfudg=1.0_jprb
    zlwfudg=1.0_jprb
  ENDIF    
ENDIF    
!print 9108,LINHOM,NHOWINH,ZSWFUDG
9108 format(1x,'LINHOM=',l8,' NHOWINH=',i2,' FUDG=',f4.2)

DO jl = kidia,kfdia
  zfcup(jl,klev+1) = _zero_
  zfcdwn(jl,klev+1) = replog
  zfsup(jl,klev+1) = _zero_
  zfsdwn(jl,klev+1) = replog
  zflux(jl,1,klev+1) = _zero_
  zflux(jl,2,klev+1) = _zero_
  zfluc(jl,1,klev+1) = _zero_
  zfluc(jl,2,klev+1) = _zero_
  zfsdnn(jl) = _zero_
  zfsdnv(jl) = _zero_
  zfcdnn(jl) = _zero_
  zfcdnv(jl) = _zero_
  zfsupn(jl) = _zero_
  zfsupv(jl) = _zero_
  zfcupn(jl) = _zero_
  zfcupv(jl) = _zero_
  zpsol(jl) = paph(jl,klev+1)
  zpmb(jl,1) = zpsol(jl) / 100._jprb
  zdt0(jl) = pts(jl) - pth(jl,klev+1)
  psudu(jl) = _zero_
  puvdf(jl) = _zero_
  pparf(jl) = _zero_
  zsudu(jl) = _zero_
  ibas(jl) = int( 0.01_jprb + pnbas(jl) )
  itop(jl) = int( 0.01_jprb + pntop(jl) )
ENDDO

DO jk=1,klev+1
  DO jl=kidia,kfdia
    cpfls(jl,jk)  = _zero_
    cpflt(jl,jk)  = _zero_
    cpfcs(jl,jk)  = _zero_
    cpfct(jl,jk)  = _zero_
    cpfdct(jl,jk) = _zero_
    cpfuct(jl,jk) = _zero_
    cpfdlt(jl,jk) = _zero_
    cpfult(jl,jk) = _zero_
    cpfdcs(jl,jk) = _zero_
    cpfucs(jl,jk) = _zero_
    cpfdls(jl,jk) = _zero_
    cpfuls(jl,jk) = _zero_
  ENDDO
ENDDO

DO jl = kidia,kfdia
  cpfrsod(jl) = _zero_
  cpemit(jl) = _zero_ 
  cpsudu(jl) = _zero_
  cpuvdf(jl) = _zero_
  cpparf(jl) = _zero_
END DO   


!*         1.1    INITIALIZE VARIOUS FIELDS
!                 -------------------------


DO jsw=1,nsw
  DO jl = kidia,kfdia
    zalbd(jl,jsw)=palbd(jl,jsw)
    zalbp(jl,jsw)=palbp(jl,jsw)
  ENDDO
ENDDO
DO jl = kidia,kfdia
  zemis(jl)  =pemis(jl)
  zemiw(jl)  =pemiw(jl)
  zmu0(jl)   =pmu0(jl)
  zuvdf(jl)  = _zero_
  zsudu(jl)  = _zero_
  zparf(jl)  = _zero_
ENDDO

DO jk = 1 , klev
  jkp1 = jk + 1
  jkl = klev+ 1 - jk
  jklp1 = jkl + 1
  DO jl = kidia,kfdia
    zpmb(jl,jk+1)=paph(jl,jkl)/100._jprb
    zoz(jl,jk)   = pozon(jl,jkl) * 46.6968_jprb / rg
    zozon(jl,jk) = pozon(jl,jkl)
    zcld0(jl,jk) = _zero_
    zfcup(jl,jk) = _zero_
    zfcdwn(jl,jk) = _zero_
    zfsup(jl,jk) = _zero_
    zfsdwn(jl,jk) = _zero_
    zflux(jl,1,jk) = _zero_
    zflux(jl,2,jk) = _zero_
    zfluc(jl,1,jk) = _zero_
    zfluc(jl,2,jk) = _zero_
  ENDDO
ENDDO


!** INPUTS ARE FULL LEVEL TEMPERATURES + SURFACE TEMPERATURE
!        INTERPOLATION TO GET HALF-LEVEL TEMPERATURES FOLLOWS
!        WHAT IS DONE IN *RADINT* AND *RADHEAT*

!* LLINTRP=.T.  Half-level temperatures on the coarse grid are
!               vertically interpolated linearly with horizontal
!               sampled pressure from the full-level temperatures
!               of the sampled grid.

!* LLINTRP=.F.  Half-level temperatures are those horizontally
!               sampled on the coarse grid

llintrp=.false.
IF (llintrp) THEN
  DO jk=2,klev
    DO jl=kidia,kfdia
      pth(jl,jk)=(pt(jl,jk-1)*pap(jl,jk-1)&
       &*(pap(jl,jk)-paph(jl,jk))&
       &+pt(jl,jk)*pap(jl,jk)*(paph(jl,jk)-pap(jl,jk-1)))&
       &*(_one_/(paph(jl,jk)*(pap(jl,jk)-pap(jl,jk-1))))
    ENDDO
  ENDDO
  IF (ltempds) THEN
    DO jl=kidia,kfdia
      pth(jl,1)= pt(jl,1)-pap(jl,1)*(pt(jl,1)-pth(jl,2))&
        &/(pap(jl,1)-paph(jl,2)) 
      pth(jl,klev+1)=pt(jl,klev)&
        &            +(paph(jl,klev+1)-pap(jl,klev))&
        &            *(pt(jl,klev)-pth(jl,klev))&
        &            /(pap(jl,klev)-paph(jl,klev))
    ENDDO
  ELSE      
    DO jl=kidia,kfdia
      pth(jl,1)= pt(jl,1)-pap(jl,1)*(pt(jl,1)-pth(jl,2))&
        &/(pap(jl,1)-paph(jl,2)) 
      pth(jl,klev+1)= pts(jl)
    ENDDO
  ENDIF    
ENDIF

DO jk=1,klev
  jkl=klev+1-jk
  jklp1=jkl+1
  DO jl=kidia,kfdia
    ztl(jl,jk)=pth(jl,jklp1)
    ztave(jl,jk)=pt(jl,jkl)
  ENDDO
ENDDO
DO jl=kidia,kfdia
  ztl(jl,klev+1)= pth(jl,1)
  zpmb(jl,klev+1) = paph(jl,1)/100._jprb
ENDDO
!***

!     ------------------------------------------------------------------

!*         2.     CLOUD AND AEROSOL PARAMETERS
!                 ----------------------------

nboxl=1
IF (kbox.EQ.1) THEN
  CALL col2box &
   & ( kidia, kfdia, klon, klev, nbox, novlp &
   & , pfrcl, pclbx &
   & )
  nboxl=nbox
END IF
zwght=1./float(nboxl) 
       
!-- initialise box-type outputs OLR, ASW, SDLW, SDSW, TAU
DO icbox=1,nboxl
  DO jl=kidia,kfdia
    olrbox(jl,icbox)=_zero_             
    aswbox(jl,icbox)=_zero_             
    slwbox(jl,icbox)=_zero_             
    sswbox(jl,icbox)=_zero_ 
    taubox(jl,icbox)=_zero_ 
  END DO  
END DO             

DO icbox=1,nboxl
  IF (kbox.EQ.1) THEN
    DO jk=1,klev
      DO jl=kidia,kfdia
        pclfr(jl,jk)=pclbx(jl,icbox,jk)
      END DO
    END DO
    
  ELSE       
    DO jk=1,klev
      DO jl=kidia,kfdia
        pclfr(jl,jk)=pfrcl(jl,jk)
      END DO
    END DO
  END IF  
  DO jl=kidia,kfdia
    psudu(jl) = _zero_
    ztauint(jl) = _zero_
  END DO  
  
!-- compute total cloud cover for that particular calculation
  DO jl=kidia,kfdia
    zcol(jl)=1.-pclfr(jl,1)
  END DO
  DO jk=2,klev
    DO jl=kidia,kfdia  
      zcol(jl)=zcol(jl)*(1.-max(pclfr(jl,jk),pclfr(jl,jk-1))) &
       & /(1.-min(pclfr(jl,jk-1),1.-repclc))
    END DO
  END DO
  DO jl=kidia,kfdia
    zcol(jl)=1.-zcol(jl)
  END DO
  





DO jk = 1 , klev
  ikl = klev + 1 - jk

!          2.1    INITIALIZE OPTICAL PROPERTIES TO CLEAR SKY VALUES
!                 -------------------------------------------------

  DO jsw = 1,nsw
    DO jl = kidia,kfdia
      ztau(jl,jsw,jk)  = _zero_
      zomega(jl,jsw,jk)= _one_
      zcg(jl,jsw,jk)   = _zero_
    ENDDO
  ENDDO
  DO jl = kidia,kfdia
    zcldsw(jl,jk)  = _zero_
    zcldld(jl,jk)  = _zero_
    zcldlu(jl,jk)  = _zero_
  ENDDO


!          2.2    CLOUD ICE AND LIQUID CONTENT AND PATH
!                 -------------------------------------

  DO jl = kidia,kfdia
!    PCLFR(JL,IKL)=MAX(REPSC,MIN(PCLFR(JL,IKL),_ONE_-REPSC))
    pclfr(jl,ikl)=max( _zero_ ,min( pclfr(jl,ikl), _one_ ))

! --- LIQUID WATER CONTENT (g.m-3) AND LIQUID WATER PATH (g.m-2)
    zlwgkg=max(pqlwp(jl,ikl)*1000._jprb,_zero_)
    ziwgkg=max(pqiwp(jl,ikl)*1000._jprb,_zero_)
!!    IF (PCLFR(JL,IKL) > (_TWO_*REPCLC)) THEN
!!      ZLWGKG=ZLWGKG/PCLFR(JL,IKL)
!!      ZIWGKG=ZIWGKG/PCLFR(JL,IKL)
!!  IF (PCLFR(JL,IKL) > REPCLC) THEN
    IF (pclfr(jl,ikl) > 15.e-06_jprb) THEN
      zlwgkg=zlwgkg/pfrcl(jl,ikl)
      ziwgkg=ziwgkg/pfrcl(jl,ikl)
    ELSE
      zlwgkg=_zero_
      ziwgkg=_zero_
    ENDIF

! --- RAIN LIQUID WATER CONTENT (g.m-3) AND LIQUID WATER PATH (g.m-2)
!    IF (PRAINT(JL,IKL).GT.(2.*REPCLC)) THEN
!      ZRWGKG=MAX(PQRAIN(JL,IKL)*1000., 0.0)
!      ZRAINT(JL)=PRAINT(JL,IKL)*3600.*1000.
!- no radiative effect of rain (for the moment)
!      ZRWGKG=0.
!      ZRAINT(JL)=0.
! ===========================================================

! Modifications Martin et al.
!    ELSE
    zrwgkg=_zero_
    zraint(jl)=_zero_
!    END IF

    IF (ibas(jl) /= 1.AND. itop(jl) /= 1 ) THEN
      zcfudg(jl)=zmtconv
    ELSE
      zcfudg(jl)=zmtfudg
    ENDIF

    zdpog=pdp(jl,ikl)/rg
    zflwp(jl)= zlwgkg*zdpog
    zfiwp(jl)= ziwgkg*zdpog
    zfrwp(jl)= zrwgkg*zdpog
    zpodt=pap(jl,ikl)/(rd*pt(jl,ikl))
    zlwc(jl)=zlwgkg*zpodt
    ziwc(jl)=ziwgkg*zpodt
!    ZRWC(JL)=ZRWGKG*ZPODT

! --- EFFECTIVE RADIUS FOR WATER, ICE AND RAIN PARTICLES

    IF (nradlp.EQ.0) THEN
! very old parametrization as f(pressure)
      zradlp(jl)=10._jprb + (100000._jprb-pap(jl,ikl))*3.5e-04_jprb

    ELSE IF (nradlp.EQ.1) THEN
! old simple distinction between land (10) and ocean (13)
      IF (plsm(jl) < _half_) THEN
        zradlp(jl)=13._jprb
      ELSE
        zradlp(jl)=10._jprb
      ENDIF

    ELSE IF (nradlp.EQ.2) THEN
!--  based on Martin et al., 1994, JAS
      IF (plsm(jl) < _half_) THEN
        zasea=150._jprb
        zd=0.33_jprb
        zntot=-1.15e-03_jprb*zasea*zasea+0.963_jprb*zasea+5.30_jprb
      ELSE
        zalnd=900._jprb
!        ZALND=600._JPRB
!        ZALND=300._JPRB
!        ZALND=1200._JPRB
        zd=0.43_jprb
        zntot=-2.10e-04_jprb*zalnd*zalnd+0.568_jprb*zalnd-27.9_jprb
      ENDIF
      
      znum=3._jprb*zlwc(jl)*(1._jprb+3._jprb*zd*zd)**2
      zden=4._jprb*rpi*zntot*(1._jprb+zd*zd)**3
      zradlp(jl)=100.*(znum/zden)**0.333_jprb
      
! 9001   format(1x,I3,1E13.5,F5.0,F5.2,f8.2,3E13.5)
      zradlp(jl)=max(zradlp(jl), 4._jprb)
      zradlp(jl)=min(zradlp(jl),16._jprb)
    END IF  
!    print *,'ZRADLP(JL) for JK=',JK,ZRADLP(JL)

! ===========================================================
! ___________________________________________________________

! rain drop from          : unused as ZRAINT is 0.
!    ZRADRD(JL)=500._JPRB*ZRAINT(JL)**0.22_JPRB
!    IF (ZFLWP(JL).GT.0.) THEN
!      ZRADRD(JL)=ZRADLP(JL)+ZRADRD(JL)
!    END IF   

!- ice particle effective radius =f(T) from Liou and Ou (1994)
 
    IF (pt(jl,ikl) < rtice) THEN
      ztempc=pt(jl,ikl)-rtt
    ELSE
      ztempc=rtice-rtt
    ENDIF
    
    zradip(jl)=326.3_jprb+ztempc*(12.42_jprb + ztempc*(0.197_jprb + ztempc*&
     &0.0012_jprb))
    zdesr(jl)=2._jprb*zradip(jl)
!    print *,'ZRADIP(JL) for JK=',JK,ZRADIP(JL),ZDESR(JL)
     
    IF (nradip.EQ. 0) THEN
      zradip(jl)= 40._jprb
      zdesr(jl)=2._jprb*zradip(jl)
!      print *,'ZRADIP(JL) for JK=',JK,ZRADIP(JL),ZDESR(JL)
      
    ELSE IF (nradip.EQ. 1) THEN 
!-- old formulation based on temperature (40-130microns)    
      zradip(jl)=max(zradip(jl),40._jprb)
      zdesr(jl)=2._jprb*zradip(jl)
!      print *,'ZRADIP(JL) for JK=',JK,ZRADIP(JL),ZDESR(JL)
      
    ELSE IF (nradip.EQ. 2) THEN  
!-- formulation following Jakob, Klein modifications to ice content    
      zradip(jl)=max(zradip(jl),30._jprb)
      zradip(jl)=min(zradip(jl),60._jprb)
      zdesr(jl)=2._jprb*zradip(jl)
!      print *,'ZRADIP(JL) for JK=',JK,ZRADIP(JL),ZDESR(JL)
      
!-- new Sun and Rikus, 1999  D_ice = f(T, IWC)    
    ELSE IF (nradip.EQ. 3 .AND. ziwc(jl).GT. _zero_ ) THEN
      ztempc=pt(jl,ikl)-83.15_jprb
      ztcels=pt(jl,ikl)-rtt
      zfsr = 1.2351_jprb +0.0105_jprb * ztcels
! Sun & Rikus, 1999      
!      ZLGIWC=LOG10( REPCLC + ZIWC(JL))
!      ZAIWC=26.1571_JPRB / ( ABS(ZLGIWC) **0.5995_JPRB )
!      ZBIWC=0.6402_JPRB + 0.1810_JPRB * ZLGIWC
! Sun, 2001
      zaiwc = 45.8966_jprb * ziwc(jl)**0.2214_jprb
      zbiwc = 0.7957_jprb * ziwc(jl)**0.2535_jprb
      zdesr(jl) = zfsr * (zaiwc + zbiwc*ztempc)
      zdesr(jl) = min( max( zdesr(jl), rminice ), 350._jprb)
      zradip(jl)= 0.5 * zdesr(jl)
!      print *,'ZRADIP(JL) for JK=',JK,ZRADIP(JL),ZDESR(JL)
    END IF
    
!-- ERA-15 definition of effective radii    
    IF (klwrad.EQ.2 .AND. nsw.EQ.2) THEN
      zradip(jl)=40._jprb
      zradlp(jl)=10._jprb + (100000._jprb-pap(jl,ikl))*3.5_jprb
!      ZSWFUDG=1._JPRB
!      ZLWFUDG=1._JPRB
      lowasyf=.false.      
      loifuec=.false.
      lrrtm=.false.
      zdesr(jl)=2._jprb*zradip(jl)
!      print *,'ZRADIP(JL) for JK=',JK,ZRADIP(JL),ZDESR(JL)
    END IF  
    
  ENDDO



!          2.3    CLOUD SHORTWAVE OPTICAL PROPERTIES
!                 ----------------------------------

!   -------------------------
! --+ SW OPTICAL PARAMETERS +  Water clouds after Fouquart (1987)
!   -------------------------  Ice clouds (Ebert, Curry, 1992)

  DO jsw=1,nsw
    DO jl = kidia,kfdia
      ztol=_zero_
      zgl =_zero_
      zol =_zero_
      ztoi=_zero_
      zgi =_zero_
      zoi =_zero_
      ztor=_zero_
      zgr =_zero_
      zor =_zero_
      IF (zflwp(jl)+zfiwp(jl)+zfrwp(jl)  /= _zero_) THEN
        IF (zflwp(jl)  /=  _zero_) THEN
          IF (nliqopt.NE.0 ) THEN
!-- SW: Slingo, 1989
            ztol = zflwp(jl)*(raswca(jsw)+raswcb(jsw)/zradlp(jl))
            zgl  = raswce(jsw)+raswcf(jsw)*zradlp(jl)
            zol  = 1. - raswcc(jsw)-raswcd(jsw)*zradlp(jl)
          ELSE          
!-- SW: Fouquart, 1991
            ztol = zflwp(jl)*(ryfwca(jsw)+ryfwcb(jsw)/zradlp(jl))
            zgl  = ryfwcf(jsw)
            zol  = ryfwcc(jsw)-ryfwcd(jsw)*exp(-ryfwce(jsw)*ztol)
          ENDIF 
        ENDIF

        IF (zfiwp(jl)  /=  _zero_) THEN
          IF (niceopt.LE.1) THEN
!-- SW: Ebert-Curry          
            ztoi = zfiwp(jl)*(rebcua(jsw)+rebcub(jsw)/zradip(jl))
            zgi  = rebcue(jsw)+rebcuf(jsw)*zradip(jl)
            zoi  = _one_ - rebcuc(jsw)-rebcud(jsw)*zradip(jl)
            
          ELSE IF (niceopt.EQ.2) THEN
!-- SW: Fu-Liou, 1993
            z1radi = 0.5 / zradip(jl)
            zbetai = rflaa0(jsw)+z1radi* rflaa1(jsw)
            ztoi = zfiwp(jl) * zbetai
            zomgi= rflbb0(jsw)+zradip(jl)*(rflbb1(jsw) + zradip(jl) &
             &   *(rflbb2(jsw)+zradip(jl)* rflbb3(jsw) ))            
            zoi  = _one_ - zomgi
            zomgp= rflcc0(jsw)+zradip(jl)*(rflcc1(jsw) + zradip(jl) &
             &   *(rflcc2(jsw)+zradip(jl)* rflcc3(jsw) )) 
            zfdel= rfldd0(jsw)+zradip(jl)*(rfldd1(jsw) + zradip(jl) &
             &   *(rfldd2(jsw)+zradip(jl)* rfldd3(jsw) )) 
            zgi  = ((1.-zfdel)*zomgp + zfdel*3.) / 3.    
                    
          ELSE IF (niceopt.EQ.3) THEN
!-- SW: Fu 1996
            z1radi = _one_ / zdesr(jl)
            zbetai = rfuaa0(jsw)+z1radi* rfuaa1(jsw)
            ztoi = zfiwp(jl) * zbetai
            zomgi= rfubb0(jsw)+zdesr(jl)*(rfubb1(jsw) + zdesr(jl) &
             &   *(rfubb2(jsw)+zdesr(jl)* rfubb3(jsw) ))            
            zoi  = _one_ - zomgi
            zgi  = rfucc0(jsw)+zdesr(jl)*(rfucc1(jsw) + zdesr(jl) &
             &   *(rfucc2(jsw)+zdesr(jl)* rfucc3(jsw) )) 
             
          ENDIF
        ENDIF

!        IF (ZFRWP(JL) .NE. 0.) THEN
!          ZTOR= ZFRWP(JL)*0.003_JPRB*_JPRBZRAINT(JL)**(-0.22_JPRB)         
!          ZOR = 1._JPRB - RROMA(JSW)*ZRAINT(JL)**RROMB(JSW)
!          ZGR = RRASY(JSW)
!        END IF   

!  - MIX of WATER and ICE CLOUDS
!        ZTAUMX= ZTOL + ZTOI + ZTOR
!        ZOMGMX= ZTOL*ZOL + ZTOI*ZOI + ZTOR*ZOR
!        ZASYMX= ZTOL*ZOL*ZGL + ZTOI*ZOI*ZGI + ZTOR*ZOR*ZGR
!
!        ZASYMX= ZASYMX/ZOMGMX
!        ZOMGMX= ZOMGMX/ZTAUMX

        IF (.NOT.linhom .OR. (linhom .AND. nhowinh.EQ.1) ) THEN
          zvl=zswfudg
          zvi=zswfudg
          zvr=0.
          ztaumx= ztol*zvl + ztoi*zvi + ztor*zvr
          zomgmx= ztol*zvl*zol + ztoi*zvi*zoi + ztor*zvr*zor
          zasymx= ztol*zvl*zol*zgl + ztoi*zvi*zoi*zgi + ztor*zvr*zor*zgr
          zasymx= zasymx/zomgmx
          zomgmx= zomgmx/ztaumx
        ELSE IF (linhom .AND. nhowinh.EQ.2) THEN
          zvl=psqlw(jl,ikl)
          zvi=psqiw(jl,ikl)
          zvr=0.
          ztaumx= ztol*zvl + ztoi*zvi + ztor*zvr
          zomgmx= ztol*zvl*zol + ztoi*zvi*zoi + ztor*zvr*zor
          zasymx= ztol*zvl*zol*zgl + ztoi*zvi*zoi*zgi + ztor*zvr*zor*zgr
          zasymx= zasymx/zomgmx
          zomgmx= zomgmx/ztaumx
        ELSE IF (linhom .AND. nhowinh.EQ.3) THEN
          zvl=prlvrl(jl,ikl)
          zvi=prlvri(jl,ikl)
          zvr=0.
          ztol1 = ztol/(1.+zvl)
          ztoi1 = ztoi/(1.+zvi)
          ztor1 = ztor/(1.+zvr) 
          ztaumx= ztol1 + ztoi1 + ztor1
          zoi=zoi/(1.+zvi*(1.-zoi))
          zgi=zgi*(1.+zvi*(1.-zoi))/(1.+zvi*(1.-zoi*zgi))
          zol=zol/(1.+zvl*(1.-zol))
          zgl=zgl*(1.+zvl*(1.-zol))/(1.+zvl*(1.-zol*zgl))
          
          zomgmx= ztol1*zol + ztoi1*zoi + ztor1*zor
          zasymx= ztol1*zol*zgl + ztoi1*zoi*zgi + ztor1*zor*zgr
          zasymx= zasymx/zomgmx
          zomgmx= zomgmx/ztaumx
        END IF  
!        print 9009,JK,JL,JSW,ZSWFUDG,PSQLW(JL,IKL),PSQIW(JL,IKL) &
!         & , PRLVRL(JL,IKL),PRLVRI(JL,IKL),ZTOL,ZOL,ZGL,ZTOI,ZOI,ZGI &
!         & , ZTAUMX,ZOMGMX,ZASYMX
9009    format(1x,3i3,14e13.6)         
        
! --- SW FINAL CLOUD OPTICAL PARAMETERS

        zcldsw(jl,jk)  = pclfr(jl,ikl)
        ztau(jl,jsw,jk)  = ztaumx
        zomega(jl,jsw,jk)= zomgmx
        zcg(jl,jsw,jk)   = zasymx
      ENDIF

! #DB             jkjlsw = 0
! #DB   IF (ZTAU(JL,JSW,JK)    .LT.00..OR.ZTAU(JL,JSW,JK)    .GT.75. .OR. &
! #DB &     (k2iii(JL)         .EQ.kio.AND.k2jjj(JL)         .EQ.kjo)     ) THEN
! #DB     IF (mod(jkjlsw,20).EQ.0)                                               &
! #DB &     write(6,575)      NLIQOPT,NICEOPT
! #DB               575  format('IN   RADLSW: CLOUD SHrtWAVE OPTICAL PROPERTIES                 '  &
! #DB &                        ,3x,'   NLIQOPT =',I3,'   NICEOPT =',I3,/                           &
! #DB &                        ,'    i    j   JL   JK',7x,'ZTAU',5x,'ZCLDSW',6x,'ZDESR'            &
! #DB &                        ,5x,'PRLVRL',5x,'PRLVRI',6x,'PQIWP',6x,'PQLWP',3x,'JSW')
! #DB             jkjlsw=jkjlsw+1
! #DB       write(6,603) k2iii(JL),k2jjj(JL),JL,IKL,ZTAU(JL,JSW,JK) ,ZCLDSW(JL,JK) , ZDESR(JL)     &
! #DB &                                                             ,PRLVRL(JL,IKL),PRLVRI(JL,IKL) &
! #DB &                                                              ,PQIWP(JL,IKL), PQLWP(JL,IKL),JSW
! #DB               603  format(4i5,7e11.3,I6)
! #DB   ENDIF

    ENDDO
  ENDDO
  
  DO jl=kidia,kfdia
    ztauint(jl)=ztauint(jl)+ztau(jl,1,jk)
  END DO  
  
  
!JL=KIDIA
!print 9109,JK,ZCLDSW(JL,JK),ZRADLP(JL),ZRADIP(JL) &
!  & , (ZTAU(JL,JSW,JK),ZOMEGA(JL,JSW,JK),ZCG(JL,JSW,JK),JSW=1,NSW)
9109 format(1x,'ClOptProp: ',i2,f7.4,2f6.1,6(1x,f7.2,1x,f7.4,1x,f6.3))  
!print *,'Radlsw after SW cloud optical properties for level JK=',JK



!          2.4    CLOUD LONGWAVE OPTICAL PROPERTIES FOR EC-OPE
!                 --------------------------------------------

!   -------------------------
! --+ LW OPTICAL PARAMETERS +  Water (and Ice) from Smith and Shi (1992)
!   -------------------------  Ice clouds (Ebert, Curry, 1992)

  IF (.NOT.lrrtm) THEN

    DO jl = kidia,kfdia
      zalfice(jl)=_zero_
      zgamice(jl)=_zero_
      zbice(jl)=_zero_
      ztice(jl)=(pt(jl,ikl)-tstand)/tstand
      zbicfu(jl)=_zero_
      zkicfu1(jl)=_zero_
      zkicfu2(jl)=_zero_
    ENDDO
    
    DO jnu= 1,nsil
      DO jl = kidia,kfdia
        zres(jl)  = xp(1,jnu)+ztice(jl)*(xp(2,jnu)+ztice(jl)*(xp(3,&
         &jnu)&
         &+ztice(jl)*(xp(4,jnu)+ztice(jl)*(xp(5,jnu)+ztice(jl)*(xp(6,&
         &jnu)&
         &)))))
        zbice(jl) = zbice(jl) + zres(jl)
        zgamice(jl) = zgamice(jl) + rebcui(jnu)*zres(jl)
        zalfice(jl) = zalfice(jl) + rebcuj(jnu)*zres(jl)
      ENDDO
    ENDDO
        
!-- Fu et al. (1998) with M'91 LW scheme    
    DO jrtm=1,16
      DO jl=kidia,kfdia
        IF (pt(jl,ikl) < 339._jprb .AND. pt(jl,ikl) >= 160._jprb) THEN
          indlay=pt(jl,ikl)-159._jprb
          ztblay =pt(jl,ikl)-int(pt(jl,ikl))
        ELSE IF (pt(jl,ikl) >= 339._jprb ) THEN
          indlay=180
          ztblay =pt(jl,ikl)-339._jprb
        ELSE IF (pt(jl,ikl) < 160._jprb) THEN
          indlay=1
          ztblay =pt(jl,ikl)-160._jprb
        END IF      
        zaddplk = totplnk(indlay+1,jrtm)-totplnk(indlay,jrtm)
        zplanck = delwave(jrtm) * (totplnk(indlay,jrtm) + ztblay*zaddplk)
        zbicfu(jl) = zbicfu(jl) + zplanck
        
        IF (ziwc(jl) > _zero_ ) THEN
! ice cloud spectral emissivity a la Fu & Liou (1993)
          zratio= 0.5 / zradip(jl)
          zmsaid = rfulio(jrtm,1) + zratio&
             &*(rfulio(jrtm,2) + zratio*rfulio(jrtm,3))
          zkicfu1(jl) = zkicfu1(jl)+ zmsaid*zplanck
          
! ice cloud spectral emissivity a la Fu et al (1998)
          z1radi = _one_ / zdesr(jl)
          zmsaid = rfueta(jrtm,1) + z1radi&
             &*(rfueta(jrtm,2) + z1radi*rfueta(jrtm,3))
          zkicfu2(jl) = zkicfu2(jl)+ zmsaid*zplanck
        END IF  
      END DO
    END DO
            
    DO jl = kidia,kfdia
      zgamice(jl) = zgamice(jl) / zbice(jl)
      zalfice(jl) = zalfice(jl) / zbice(jl)
      zkicfu1(jl) = zkicfu1(jl) / zbicfu(jl)
      zkicfu2(jl) = zkicfu2(jl) / zbicfu(jl)
      
      IF (zflwp(jl)+zfiwp(jl) /= _zero_) THEN

        IF (klwrad.EQ.2) THEN        
! ice cloud emissivity a la Smith-Shi
          zmulti=1.2_jprb-0.006_jprb*zradip(jl)
          zmsaid= 0.113_jprb*zmulti
          zmsaiu= 0.093_jprb*zmulti
          zmultl=1.2_jprb-0.006_jprb*zradlp(jl)
          zmsald= 0.158_jprb*zmultl
          zmsalu= 0.130_jprb*zmultl
          zzflwp= zflwp(jl)
          zzfiwp= zfiwp(jl)
          
        ELSE IF (klwrad.EQ.0) THEN  
          
          IF (nliqopt.EQ.0) THEN
! water cloud emissivity a la Smith & Shi (1992)
            zmultl=1.2_jprb-0.006_jprb*zradlp(jl)
            zmsald= 0.158_jprb*zmultl
            zmsalu= 0.130_jprb*zmultl
          
          ELSE
! water cloud emissivity a la Savijarvi (1997)
            zmsalu= 0.2441_jprb-0.0105_jprb*zradlp(jl)
            zmsald= 1.2154_jprb*zmsalu
          
          END IF  
          
          IF (niceopt.EQ.0) THEN          
! ice cloud emissivity a la Smith & Shi (1992)
            zmulti=1.2_jprb-0.006_jprb*zradip(jl)
            zmsaid= 0.113_jprb*zmulti
            zmsaiu= 0.093_jprb*zmulti

          ELSE IF (niceopt.EQ.1) THEN
! ice cloud emissivity a la Ebert & Curry (1992)
            zmsaid= 1.66_jprb*(zalfice(jl)+zgamice(jl)/zradip(jl))
            zmsaiu= zmsaid
         
          ELSE IF (niceopt.EQ.2) THEN  
! ice cloud emissivity a la Fu & Liou (1993)
            zmsaid= 1.66_jprb*zkicfu1(jl)
            zmsaiu= zmsaid
          
          ELSE IF (niceopt.EQ.3) THEN  
! ice cloud emissivity a la Fu et al. (1998)
            zmsaid= 1.66_jprb*zkicfu2(jl)
            zmsaiu= zmsaid
          END IF  
          
! introduce inhomogeneity factor also in LW         
          zzflwp= zflwp(jl) * zlwfudg
          zzfiwp= zfiwp(jl) * zlwfudg
        END IF
          
! effective cloudiness accounting for condensed water
        zcldld(jl,jk) = pclfr(jl,ikl)*(_one_-exp(-zmsald*zzflwp-zmsaid* &
          &zzfiwp))
        zcldlu(jl,jk) = pclfr(jl,ikl)*(_one_-exp(-zmsalu*zzflwp-zmsaiu* &
          &zzfiwp))
          
      END IF    
    ENDDO
    
!  print *,'Radlsw after LW0 cloud optical properties for level JK=',JK

  ELSE

!          2.5    CLOUD LONGWAVE OPTICAL PROPERTIES FOR RRTM
!                 ------------------------------------------

!   -------------------------
! --+ LW OPTICAL PARAMETERS +  Water (and Ice) from Savijarvi (1998)
!   -------------------------  Ice clouds (Ebert, Curry, 1992)

! No need for a fixed diffusivity factor, accounted for spectrally below
! The detailed spectral structure does not require defining upward and
! downward effective optical properties

! #DB             jkjllw=0

    DO jrtm=1,16
      DO jl = kidia,kfdia
        ztaucld(jl,jk,jrtm) = _zero_
        zmsald = _zero_
        zmsaid = _zero_
        
        IF (zflwp(jl)+zfiwp(jl) /= _zero_) THEN
    
          IF (nliqopt.EQ.0) THEN
! water cloud total emissivity a la Smith and Shi (1992)
            zmultl=1.2_jprb-0.006_jprb*zradlp(jl)
            zmsald= 0.144_jprb*zmultl / 1.66_jprb
            
          ELSE IF (nliqopt.EQ.1) THEN
! water cloud spectral emissivity a la Savijarvi (1997)
            zmsald= rhsavi(jrtm,1) + zradlp(jl)&
             &*(rhsavi(jrtm,2) + zradlp(jl)*rhsavi(jrtm,3))
             
          ELSE IF (nliqopt.EQ.2) THEN
! water cloud spectral emissivity a la Lindner and Li (2000)
            z1radl = _one_ / zradlp(jl)
!            ZMSALD = RLINLI(JRTM,1) + Z1RADL*(RLINLI(JRTM,2) + Z1RADL*&
!            &       (RLINLI(JRTM,3) + Z1RADL*(RLINLI(JRTM,4) + Z1RADL*&
!            &        RLINLI(JRTM,5) )))
            
            zmsald = rlinli(jrtm,1)+zradlp(jl)*rlinli(jrtm,2)+ z1radl*&
            &       (rlinli(jrtm,3) + z1radl*(rlinli(jrtm,4) + z1radl*&
            &        rlinli(jrtm,5) ))
          
          END IF  

          IF (niceopt.EQ.0) THEN
! ice cloud emissivity a la Smith & Shi (1992)
            zmulti=1.2_jprb-0.006_jprb*zradip(jl)
            zmsaid= 0.108_jprb*zmulti / 1.66_jprb
                   
          ELSE IF (niceopt.EQ.1) THEN
! ice cloud spectral emissivity a la Ebert-Curry (1992)
            zmsaid= rebcuh(jrtm)+rebcug(jrtm)/zradip(jl)
            
          ELSE IF (niceopt.EQ.2) THEN
! ice cloud spectral emissivity a la Fu & Liou (1993)
            zratio= 0.5 / zradip(jl)
            zmsaid = rfulio(jrtm,1) + zratio&
             &*(rfulio(jrtm,2) + zratio*rfulio(jrtm,3))
             
          ELSE IF (niceopt.EQ.3) THEN
! ice cloud spectral emissivity a la Fu et al (1998)
            z1radi = _one_ / zdesr(jl)
            zmsaid = rfueta(jrtm,1) + z1radi&
             &*(rfueta(jrtm,2) + z1radi*rfueta(jrtm,3))
             
          END IF    

          IF (.NOT.linhom .OR. (linhom .AND. nhowinh.EQ.1) ) THEN
            zvl=zlwfudg
            zvi=zlwfudg
          ELSE IF (linhom .AND. nhowinh.EQ.2) THEN
            zvl=psqlw(jl,ikl)
            zvi=psqiw(jl,ikl)
          ELSE IF (linhom .AND. nhowinh.EQ.3) THEN
            zvl=_one_/(_one_+prlvrl(jl,ikl))
            zvi=_one_/(_one_+prlvri(jl,ikl))
          END IF  
          
          ztaud = zvl*zmsald*zflwp(jl)+zvi*zmsaid*zfiwp(jl)

! #DB     write(30,333) ZTAUD,ZVL,ZMSALD,ZFLWP(JL),ZVI,ZMSAID,ZFIWP(JL),PQIWP(JL,IKL),PQLWP(JL,IKL)
! #DB              333  format(9e14.6)

! Diffusivity correction within clouds a la Savijarvi
!          ZDIFFD=MIN(MAX(1.517_JPRB-0.156_JPRB*LOG(ZTAUD) , _ONE_) , _TWO_)

          zdiffd=1.66_jprb
          ztaucld(jl,jk,jrtm) = max(_zero_,ztaud*zdiffd)
        ENDIF
        
! #DB   IF (ZTAUCLD(JL,JK,JRTM).LT.00..OR.ZTAUCLD(JL,JK,JRTM).GT.75. .OR. &
! #DB &     (k2iii(JL)         .EQ.kio.AND.k2jjj(JL)         .EQ.kjo)     ) THEN
! #DB     IF (mod(jkjllw,20).EQ.0)                                                 &
! #DB &     write(6,600) JRTM,NLIQOPT,NICEOPT
! #DB               600  format('IN   RADLSW: CLOUD LONGWAVE OPTICAL PROPERTIES FOR RRTM, JRTM =',I3 &
! #DB &                        ,'   NLIQOPT =',I3,'   NICEOPT =',I3,/                                &
! #DB &                        ,'    i    j   JL   JK',7x,'ZTAU'  ,6x,'ZFLWP' ,6x,'ZFIWP'            &
! #DB &                        ,                       5x,'ZRADLP',5x,'ZRADIP',6x,'PQIWP',6x,'PQLWP')
! #DB             jkjllw=jkjllw+1
! #DB       write(6,601) k2iii(JL),k2jjj(JL),JL,IKL,ZTAUCLD(JL,JK,JRTM), ZFLWP(JL)    , ZFIWP(JL)    &
! #DB &                                                                ,ZRADLP(JL)    ,ZRADIP(JL)    &
! #DB &                                                                 ,PQIWP(JL,IKL), PQLWP(JL,IKL)
! #DB               601  format(4i5,10e11.3)
! #DB   ENDIF

      ENDDO
    ENDDO
!  print *,'Radlsw after LW1 cloud optical properties for level JK=',JK

  ENDIF

ENDDO

nuaer = nua
ntraer = ntra

!     ------------------------------------------------------------------

!*         2.6    DIFFUSIVITY FACTOR OR SATELLITE VIEWING ANGLE
!                 ---------------------------------------------


DO jl = kidia,kfdia
  zview(jl) = diff
  zemit(jl) = _zero_
ENDDO

!     ------------------------------------------------------------------

!*         3.     CALL LONGWAVE RADIATION CODE
!                 ----------------------------


!*         3.1    FULL LONGWAVE RADIATION COMPUTATIONS
!                 ------------------------------------

!print *,'Just before calling the radiation schemes'
!JL=KIDIA
!DO JK=1,KLEV
!  IKL=KLEV+1-JK
!  PRINT 9311,JK,PCLFR(JL,IKL),ZCLDLD(JL,JK),ZTAUCLD(JL,JK,1) &
!  & ,(ZTAU(JL,JSW,JK),ZOMEGA(JL,JSW,JK),ZCG(JL,JSW,JK),JSW=1,NSW) &
!  & ,(PAER(JL,JAE,JK),JAE=1,6)
9311 format(1x,i2,2f8.5,26e12.5)
!END DO
!print *,'KLWRAD=',KLWRAD,'  LPHYLIN: ',LPHYLIN,'  LRRTM: ',LRRTM

IF (.NOT.lphylin) THEN
  IF ( .NOT. lrrtm) THEN

      
    IF (klwrad .EQ. 2) THEN
      CALL olw &
       & ( kidia, kfdia , klon  , klev &
       & , pcco2, zcldld, zcldlu &
       & , pdp  , zdt0  , zemis  &
       & , paph , pozon , pth &
       & , paer , pt    , zview , pq &
       & , zcoolr,zcoolc, zflux, zfluc &
       & )
       
    ELSE IF (klwrad .EQ. 0) THEN  
     
      CALL lw &
       &( kidia , kfdia , klon  , klev , kmode &
       &, pcco2 , zcldld, zcldlu &
       &, pdp   , zdt0  , zemis , zemiw &
       &, zpmb  , pozon , ztl &
       &, paer  , ztave , zview , pq &
       &, zcoolr, zcoolc, zemit , zflux, zfluc &
       &)
       
     END IF  

  ELSE


!*         3.2    FULL LONGWAVE RADIATION COMPUTATIONS - RRTM
!                 ------------------------------------   ----

!  i)  pass POZN (ozone mmr concentration) to RRTM; remove pressure
!      weighting applied to POZON in driverMC (below)
!  ii) pass ZEMIS and ZEMIW to RRTM; return ZEMIT from RRTM
!  iii)pass ZTAUCLD, cloud optical depths (water+ice) to RRTM, 
!      computed from equations above
!  iv) pass ECRT arrays to RRTM arrays in interface routine ECRTATM
!      in module rrtm_ecrt.f

    DO jl = kidia,kfdia
      DO jk = 1, klev
        zozn(jl,jk) = pozon(jl,jk)/pdp(jl,jk)
      ENDDO
    ENDDO

! #DB             jkjllw = 0
! #DB DO JL = KIDIA,KFDIA
! #DB DO JK = 1, KLEV
! #DB   IKL = KLEV + 1 - JK
! #DB      JAERmin=1
! #DB      JAERmax=6
! #DB      PAERmin=1000.
! #DB      PAERmax=0.
! #DB     jTAUCLDmin=1
! #DB     jTAUCLDmax=16
! #DB      TAUCLDmin=1000.
! #DB      TAUCLDmax=0.
! #DB   DO JRTM=1,16
! #DB    IF (ZTAUCLD(JL,JK,JRTM).LT.TAUCLDmin)                      THEN
! #DB      jTAUCLDmin=JRTM
! #DB       TAUCLDmin=ZTAUCLD(JL,JK,JRTM)
! #DB    END IF
! #DB    IF (ZTAUCLD(JL,JK,JRTM).GT.TAUCLDmax)                      THEN
! #DB      jTAUCLDmax=JRTM
! #DB       TAUCLDmax=ZTAUCLD(JL,JK,JRTM)
! #DB    END IF
! #DB   ENDDO
! #DB   DO JAE =1,6
! #DB    IF (PAER(JL,JAE,JK).LT.PAERmin)                            THEN
! #DB      JAERmin=JAE
! #DB      PAERmin=PAER(JL,JAE,JK)
! #DB    END IF
! #DB    IF (PAER(JL,JAE,JK).GT.PAERmax)                            THEN
! #DB      JAERmax=JAE
! #DB      PAERmax=PAER(JL,JAE,JK)
! #DB    END IF
! #DB   ENDDO
! #DB   IF (TAUCLDmin.LT.0..OR.TAUCLDmax.GT.75.)                    THEN
! #DB     IF (mod(jkjllw,20).EQ.0)                                                         &
! #DB &     write(6,515)
! #DB               515  format('IN   RADLSW: BEFORE RRTM_RRTM_140GP CALL',/               &
! #DB &                         ,'    i    j   JL   JK',7x,'ZOZN',5x,'ZCLDSW'              &
! #DB &                         ,4x,'ZTAUCLDmin',4x,'ZTAUCLDmax'                           &
! #DB &                         ,4x,'PAERmin',4x,'PAERmax',6x,'PQIWP',6x,'PQLWP',9x,'PQ')
! #DB             jkjllw = jkjllw + 1
! #DB       write(6,602) k2iii(JL),k2jjj(JL),JL,JK,ZOZN(JL,JK),ZCLDSW(JL,JK)               &
! #DB &                                           ,jTAUCLDmin,ZTAUCLD(JL,JK,jTAUCLDmin)    &
! #DB &                                           ,jTAUCLDmax,ZTAUCLD(JL,JK,jTAUCLDmax)    &
! #DB &                                           ,PAER(JL,JAERmin,JK),PAER(JL,JAERmax,JK) &
! #DB &                                           ,PQIWP(JL,IKL),PQLWP(JL,IKL),PQ(JL,IKL)
! #DB               602  format(4i5,2e11.3,2(i3,e11.3),8e11.3)
! #DB   ENDIF
! #DB ENDDO
! #DB ENDDO

!    print *,'Just before calling RRTM'

    CALL rrtm_rrtm_140gp &
     &( kidia , kfdia , klon  , klev &
     &, paer  , paph  , pap   &
     &, pts   , pth   , pt     &
     &, zemis , zemiw &
     &, pq    , pcco2 , zozn  , zcldsw  , ztaucld &
     &, zemit , zflux , zfluc , ztclear &
     &)
     
!     print *,'just after RRTM'    

  ENDIF
ELSE
  zcoolr(:,:) = _zero_
  zcoolc(:,:) = _zero_
  zemit(:)   = _zero_
  zflux(:,:,:)= _zero_
  zfluc(:,:,:)= _zero_
ENDIF

!     ------------------------------------------------------------------

!*         4.     CALL SHORTWAVE RADIATION CODE
!                 -----------------------------


zrmuz=_zero_
DO jl = kidia,kfdia
  zrmuz = max(zrmuz, zmu0(jl))
ENDDO

IF (zrmuz > _zero_) THEN
!print *,'CALL SW'    

  CALL sw &
   &( kidia , kfdia , klon  , klev  , kaer &
   &, prii0 , pcco2 , zpsol , zalbd , zalbp , pq   , pqs &
   &, zmu0  , zcg   , zcldsw, pdp   , zomega, zoz  , zpmb &
   &, ztau  , ztave , paer &
   &, zheatr, zfsdwn, zfsup , zheatc, zfcdwn, zfcup &
   &, zfsdnn, zfsdnv, zfsupn, zfsupv &
   &, zfcdnn, zfcdnv, zfcupn, zfcupv &
   &, zsudu , zuvdf , zparf &
   &)
   
!     print *,'just after SW'    
!     JL=KIDIA
!     print *,'just after SW UV & PAR ',ZUVDF(JL),ZPARF(JL)

ENDIF

! #DB             jkjlsw = 0
! #DB DO JL = KIDIA,KFDIA
! #DB DO JK = 1,KLEV
! #DB   IF (k2iii(JL).EQ.kio.AND.k2jjj(JL).EQ.kjo) THEN
! #DB     IF (mod(jkjlsw,20).EQ.0)                                                &
! #DB       write(6,525)
! #DB               525 format('IN   RADLSW: AFTER  SW              CALL',/       &
! #DB &                       ,'    i    j   JL   JK'                             &
! #DB &                       ,4x,'ZFCDWN',5x,'ZFCUP' ,4x,'ZFSDNN',4x,'ZFCDNN'    &
! #DB &                       ,4x,'ZFSDNV',4x,'ZFSUPN',4x,'ZFSUPV',4x,'ZFCDNN'    &
! #DB &                       ,4x,'ZFCDNV',4x,'ZFCUPN',4x,'ZFCUPV')
! #DB             jkjlsw = jkjlsw + 1
! #DB       write(6,605) k2iii(JL),k2jjj(JL),JL,JK,ZFCDWN(JL,JK),ZFCUP(JL,JK)     &
! #DB &               ,ZFSDNN(JL),ZFCDNN(JL),ZFSDNV(JL),ZFSUPN(JL),ZFSUPV(JL)     &
! #DB &               ,ZFCDNN(JL),ZFCDNV(JL),ZFCUPN(JL),ZFCUPV(JL)
! #DB               605  format(4i5,11e10.3)
! #DB   ENDIF
! #DB ENDDO
! #DB ENDDO

!     ------------------------------------------------------------------

!*         5.     FILL UP THE MODEL NET LW AND SW RADIATIVE FLUXES
!                 ------------------------------------------------


DO jkl = 1 , klev+1
  jk = klev+1 + 1 - jkl
  DO jl = kidia,kfdia
!    print 9506,JK,ZFSDWN(JL,JK),ZFSUP(JL,JK),ZFLUX(JL,1,JK),ZFLUX(JL,2,JK) &
!    & , ZFCDWN(JL,JK),ZFCUP(JL,JK),ZFLUC(JL,1,JK),ZFLUC(JL,2,JK)
9506 format(1x,i3,8f10.3)    
  
    cpfls(jl,jkl) =cpfls(jl,jkl) +zwght*(zfsdwn(jl,jk) - zfsup(jl,jk))
    cpflt(jl,jkl) =cpflt(jl,jkl) +zwght*(- zflux(jl,1,jk) - zflux(jl,2,jk))
    cpfcs(jl,jkl) =cpfcs(jl,jkl) +zwght*(zfcdwn(jl,jk)  - zfcup(jl,jk))
    cpfct(jl,jkl) =cpfct(jl,jkl) +zwght*(- zfluc(jl,1,jk) - zfluc(jl,2,jk))
    cpfdct(jl,jkl)=cpfdct(jl,jkl)+zwght*zfluc(jl,2,jk)
    cpfuct(jl,jkl)=cpfuct(jl,jkl)+zwght*zfluc(jl,1,jk)
    cpfdlt(jl,jkl)=cpfdlt(jl,jkl)+zwght*zflux(jl,2,jk)
    cpfult(jl,jkl)=cpfult(jl,jkl)+zwght*zflux(jl,1,jk)
    cpfdcs(jl,jkl)=cpfdcs(jl,jkl)+zwght*zfcdwn(jl,jk)
    cpfucs(jl,jkl)=cpfucs(jl,jkl)+zwght*zfcup(jl,jk)
    cpfdls(jl,jkl)=cpfdls(jl,jkl)+zwght*zfsdwn(jl,jk)
    cpfuls(jl,jkl)=cpfuls(jl,jkl)+zwght*zfsup(jl,jk)
  ENDDO
ENDDO

DO jl = kidia,kfdia
! print 9507,ZFSDWN(JL,1),ZSUDU(JL),ZUVDF(JL),ZPARF(JL)
9507 format(1x,'SW Global Normal UV & PAR:',5f10.3)
  
  cpfrsod(jl) = cpfrsod(jl) + zwght*zfsdwn(jl,1)
  cpemit(jl) = cpemit(jl) + zwght*zemit(jl)
  cpsudu(jl) = cpsudu(jl) + zwght*zsudu(jl)
  cpuvdf(jl) = cpuvdf(jl) + zwght*zuvdf(jl)
  cpparf(jl) = cpparf(jl) + zwght*zparf(jl)
  
  aswbox(jl,icbox) = -zfsdwn(jl,klev+1) + zfsup(jl,klev+1)
  olrbox(jl,icbox) = -zflux(jl,1,klev+1)
  slwbox(jl,icbox) = -zflux(jl,2,1)
  sswbox(jl,icbox) = -zfsdwn(jl,1)
  taubox(jl,icbox) = ztauint(jl)
  ztca(jl) = ztca(jl) + zwght*zcol(jl)
!  print 9508,ICBOX,ASWBOX(JL,ICBOX),OLRBOX(JL,ICBOX),SLWBOX(JL,ICBOX) &
!  & ,SSWBOX(JL,ICBOX),TAUBOX(JL,ICBOX),ZCOL(JL),ZTCA(JL),ZTCC(JL)
9508 format(1x,'radlsw',i3,5f10.3,1x,3f7.4)  
ENDDO


ENDDO
!
!-- end of box-type calculations
!      
 
DO jk = 1 , klev+1
  DO jl = kidia,kfdia
    pfls(jl,jk)  = cpfls(jl,jk) 
    pflt(jl,jk)  = cpflt(jl,jk) 
    pfcs(jl,jk)  = cpfcs(jl,jk) 
    pfct(jl,jk)  = cpfct(jl,jk) 
    pfdct(jl,jk) = cpfdct(jl,jk)
    pfuct(jl,jk) = cpfuct(jl,jk)
    pfdlt(jl,jk) = cpfdlt(jl,jk)
    pfult(jl,jk) = cpfult(jl,jk)
    pfdcs(jl,jk) = cpfdcs(jl,jk)
    pfucs(jl,jk) = cpfucs(jl,jk)
    pfdls(jl,jk) = cpfdls(jl,jk)
    pfuls(jl,jk) = cpfuls(jl,jk)
  ENDDO
ENDDO

DO jl = kidia,kfdia
  pfrsod(jl) = cpfrsod(jl)
  pemit(jl) = cpemit(jl)
  psudu(jl) = cpsudu(jl)
  puvdf(jl) = cpuvdf(jl)
  pparf(jl) = cpparf(jl)
ENDDO

!-- re-organize the box-tyoe output arrays in decreasing order of TAU
DO jl=kidia,kfdia
  DO icbox=2,nbox
    ztoi=taubox(jl,icbox)
    DO imov=icbox-1,1,-1
      IF(taubox(jl,imov).LE.ztoi) GO TO 8001
        taubox(jl,imov+1)=taubox(jl,imov)
    END DO
    imov=0
8001 CONTINUE
    taubox(jl,imov+1)=ztoi
  END DO  
END DO

!-- re-organize the box-type output arrays in decreasing order of ASW
DO jl=kidia,kfdia
  DO icbox=2,nbox
    zasw=aswbox(jl,icbox)
    DO imov=icbox-1,1,-1
      IF(aswbox(jl,imov).LE.zasw) GO TO 8002
        aswbox(jl,imov+1)=aswbox(jl,imov)
    END DO
    imov=0
8002 CONTINUE
    aswbox(jl,imov+1)=zasw
  END DO  
END DO

!-- re-organize the box-tyoe output arrays in decreasing order of -OLR
DO jl=kidia,kfdia
  DO icbox=2,nbox
    zolr=olrbox(jl,icbox)
    DO imov=icbox-1,1,-1
      IF(olrbox(jl,imov).LE.zolr) GO TO 8003
        olrbox(jl,imov+1)=olrbox(jl,imov)
    END DO
    imov=0
8003 CONTINUE
    olrbox(jl,imov+1)=zolr
  END DO  
END DO

!-- re-organize the box-tyoe output arrays in decreasing order of SLW
DO jl=kidia,kfdia
  DO icbox=2,nbox
    zslw=slwbox(jl,icbox)
    DO imov=icbox-1,1,-1
      IF(slwbox(jl,imov).LE.zslw) GO TO 8004
        slwbox(jl,imov+1)=slwbox(jl,imov)
    END DO
    imov=0
8004 CONTINUE
    slwbox(jl,imov+1)=zslw
  END DO  
END DO

!-- re-organize the box-type output arrays in decreasing order of -SSW
DO jl=kidia,kfdia
  DO icbox=2,nbox
    zssw=sswbox(jl,icbox)
    DO imov=icbox-1,1,-1
      IF(sswbox(jl,imov).LE.zssw) GO TO 8005
        sswbox(jl,imov+1)=sswbox(jl,imov)
    END DO
    imov=0
8005 CONTINUE
    sswbox(jl,imov+1)=zssw
  END DO  
END DO

!-- put all arrays as positive numbers for plotting
DO jl=kidia,kfdia
  DO icbox=1,nbox
    aswbox(jl,icbox)=-aswbox(jl,icbox)
    olrbox(jl,icbox)=-olrbox(jl,icbox)
    sswbox(jl,icbox)=-sswbox(jl,icbox)
  END DO  
END DO

!     --------------------------------------------------------------

RETURN
END SUBROUTINE radlsw