SUBROUTINE RADLSW &
 & ( KIDIA, KFDIA , KLON , KLEV  , KMODE, KAER,&
 & PRII0,&
 & PAER , PALBD , PALBP, PAPH , PAP,&
 & PCCNL, PCCNO,&
 & PCCO2, PCLFR , PDP  , PEMIS, PEMIW , PLSM , PMU0, POZON,&
 & PQ   , PQIWP , PQLWP, PQS  , PQRAIN, PRAINT,&
 & PTH  , PT    , PTS  , PNBAS, PNTOP,&
 & PREF_LIQ, PREF_ICE,&
 & PEMIT, PFCT  , PFLT , PFCS , PFLS,&
 & PFRSOD,PSUDU , PUVDF, PPARF, PPARCF, PTINCF,&
 & PSFSWDIR, PSFSWDIF,PFSDNN,PFSDNV ,&  
 & LRDUST,PPIZA_DST,PCGA_DST,PTAUREL_DST,&
 & PTAU_LW,&
 & PFLUX,PFLUC,PFSDN ,PFSUP , PFSCDN , PFSCUP,&
 & PTOAB, PTOACB )   !FC

use write_field_phy
! Temporary fix waiting for cleaner interface (or not)
USE clesphys_mod_h, ONLY: NSW

!**** *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
! PCCNL  : (KLON)            ; CCN CONCENTRATION OVER LAND
! PCCNO  : (KLON)            ; CCN CONCENTRATION OVER OCEAN
! PCCO2  :                   ; CONCENTRATION IN CO2 (KG/KG)
! PCLFR  : (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)       ; OZONE AMOUNT in LAYER (KG/KG*PA)
! PQ     : (KLON,KLEV)       ; SPECIFIC HUMIDITY KG/KG
! 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)
! PTH    : (KLON,KLEV+1)     ; HALF LEVEL TEMPERATURE
! PT     : (KLON,KLEV)       ; FULL LEVEL TEMPERATURE
! PTS    : (KLON)            ; SURFACE TEMPERATURE
! LDDUST                     ; Dust properties switch
! PPIZA_DST  : (KPROMA,KLEV,NSW); Single scattering albedo of dust 
! PCGA_DST   : (KPROMA,KLEV,NSW); Assymetry factor for dust 
! PTAUREL_DST: (KPROMA,KLEV,NSW); Optical depth of dust relative to at 550nm
! PTAU_LW  (KPROMA,KLEV,NLW); LW Optical depth of aerosols 
! PREF_LIQ (KPROMA,KLEV)        ; Liquid droplet radius (um)
! PREF_ICE (KPROMA,KLEV)        ; Ice crystal radius (um)
!     ==== 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
! PPARF  : (KLON)            ; PHOTOSYNTHETICALLY ACTIVE RADIATION
! PUVDF  : (KLON)            ; UV(-B) RADIATION
! PPARCF : (KLON)            ; CLEAR-SKY PHOTOSYNTHETICALLY ACTIVE RADIATION
! PTINCF : (KLON)            ; TOA INCIDENT SOLAR RADIATION 
! Ajout flux LW et SW montants et descendants, et ciel clair (MPL 19.12.08)
! PFLUX  : (KLON,2,KLEV+1)   ; LW total sky flux (1=up, 2=down)
! PFLUC  : (KLON,2,KLEV+1)   ; LW clear sky flux (1=up, 2=down)
!FC PTOAB : (klon,NLW)           ; LW total sky TOA flux bande
!FC PTOACB : (klon,NLW)          ; LW clear sky TOA flux bande
! PFSDN(KLON,KLEV+1)         ; SW total sky flux down
! PFSUP(KLON,KLEV+1)         ; SW total sky flux up
! PFSCDN(KLON,KLEV+1)        ; SW clear sky flux down
! PFSCUP(KLON,KLEV+1)        ; SW clear sky flux up



!        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
!        JJMorcrette : 010112 Sun-Rikus ice particle Diameter
!        JJMorcrette : 010301 cleaning liq/ice cloud optical properties
!        JJMorcrette : 011005 CCN --> Re liquid water clouds
!        JJMorcrette : 011108 Safety checks
!        JJMorcrette : 011108 Safety checks
!        DJSalmond   : 020211 Check before R-To-R
!        JJMorcrette : 020901 PAR & UV
!        M.Hamrud      01-Oct-2003 CY28 Cleaning
!        JJMorcrette : 050402 New sets of optical properties (NB: inactive)
!        Y.Seity       04-11-18 : add 4 arguments for AROME externalized surface
!        Y.Seity       05-10-10 : add 3 optional arg. for dust SW properties
!        JJMorcrette 20060721 PP of clear-sky PAR and TOA incident solar radiation 
!-----------------------------------------------------------------------

USE PARKIND1  ,ONLY : JPIM     ,JPRB
USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK

USE YOMCST   , ONLY : RG       ,RD       ,RTT      ,RPI
!USE YOERAD   , ONLY : NSW      ,LRRTM    ,LCCNL    ,LCCNO, LDIFFC, &
! NSW mis dans .def MPL 20140211
USE YOERAD   , ONLY : NLW, LRRTM    ,LCCNL    ,LCCNO, LDIFFC, &
 & NRADIP , NRADLP , NICEOPT, NLIQOPT, NINHOM ,NLAYINH ,&
 & RCCNLND, RCCNSEA, RLWINHF, RSWINHF, RRe2De ,&  
 & LEDBUG  
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   ,RFUETA   ,RFUETB   ,RFUETC  ,RASWCA   ,&
 & RASWCB   ,RASWCC   ,RASWCD   ,RASWCE   ,RASWCF   ,&
 & RFUAA0   ,RFUAA1   ,RFUBB0   ,RFUBB1   ,RFUBB2   ,&
 & RFUBB3   ,RFUCC0   ,RFUCC1   ,RFUCC2   ,RFUCC3   ,& 
 & RLILIA   ,RLILIB  
USE YOERDU        , ONLY : NUAER    ,NTRAER   ,REPLOG   ,REPSC    ,REPSCW   ,DIFF
!USE YOETHF        , ONLY : RTICE
USE YOEPHLI       , ONLY : LPHYLIN
USE YOERRTWN      , ONLY :                     DELWAVE   ,TOTPLNK   

USE YOMLUN_IFSAUX , ONLY : NULOUT
USE YOMCT3        , ONLY : NSTEP
!FC
USE PARRRTM, ONLY : JPGPT, GP_PER_LWBAND, JPBAND


IMPLICIT NONE

!!include "clesphys.h"
!!include "clesrrtm.h"
include "YOETHF.h"
INTEGER(KIND=JPIM),INTENT(IN)    :: KLON 
INTEGER(KIND=JPIM),INTENT(IN)    :: KLEV 
INTEGER(KIND=JPIM),INTENT(IN)    :: KIDIA 
INTEGER(KIND=JPIM),INTENT(IN)    :: KFDIA 
INTEGER(KIND=JPIM),INTENT(IN)    :: KMODE 
INTEGER(KIND=JPIM),INTENT(IN)    :: KAER 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PRII0 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PAER(KLON,6,KLEV) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PALBD(KLON,NSW) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PALBP(KLON,NSW) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PAPH(KLON,KLEV+1) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PAP(KLON,KLEV) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PCCNL(KLON) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PCCNO(KLON) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PCCO2 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PCLFR(KLON,KLEV) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PDP(KLON,KLEV) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PEMIS(KLON) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PEMIW(KLON) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PLSM(KLON) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PMU0(KLON) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: POZON(KLON,KLEV) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PQ(KLON,KLEV) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PQIWP(KLON,KLEV) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PQLWP(KLON,KLEV) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PQS(KLON,KLEV) 
REAL(KIND=JPRB)                  :: PQRAIN(KLON,KLEV) ! Argument NOT used
REAL(KIND=JPRB)                  :: PRAINT(KLON,KLEV) ! Argument NOT used
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTH(KLON,KLEV+1) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PT(KLON,KLEV) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTS(KLON) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PNBAS(KLON) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PNTOP(KLON) 
LOGICAL           ,INTENT(IN)    :: LRDUST
REAL(KIND=JPRB)   ,INTENT(IN)    :: PPIZA_DST(KLON,KLEV,NSW)
REAL(KIND=JPRB)   ,INTENT(IN)    :: PCGA_DST(KLON,KLEV,NSW)
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTAUREL_DST(KLON,KLEV,NSW)
!--C.Kleinschmitt
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTAU_LW(KLON,KLEV,NLW)
!--end
REAL(KIND=JPRB)   ,INTENT(IN)    :: PREF_LIQ(KLON,KLEV)
REAL(KIND=JPRB)   ,INTENT(IN)    :: PREF_ICE(KLON,KLEV)
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PEMIT(KLON) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFCT(KLON,KLEV+1) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFLT(KLON,KLEV+1) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFCS(KLON,KLEV+1) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFLS(KLON,KLEV+1) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFRSOD(KLON) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PSUDU(KLON) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PUVDF(KLON) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PPARF(KLON) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PPARCF(KLON), PTINCF(KLON) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PSFSWDIR(KLON,NSW) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PSFSWDIF(KLON,NSW) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFSDNN(KLON) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFSDNV(KLON) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFLUX(KLON,2,KLEV+1) ! LW total sky flux (1=up, 2=down)
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFLUC(KLON,2,KLEV+1) ! LW clear sky flux (1=up, 2=down)
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PTOAB(klon,NLW)  ! FC LW full sectral bande TOA
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PTOACB(klon,NLW) ! FC LW clear spectral TOA
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFSDN(KLON,KLEV+1)   ! SW total sky flux down
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFSUP(KLON,KLEV+1)   ! SW total sky flux up
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFSCDN(KLON,KLEV+1)  ! SW clear sky flux down
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFSCUP(KLON,KLEV+1)  ! SW clear sky flux up


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

!*       0.1   ARGUMENTS.
!              ----------
!     ==== COMPUTED IN RADLSW ===
!     -----------------------------------------------------------------

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

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

INTEGER(KIND=JPIM) :: IBAS(KLON)     , ITOP(KLON)

REAL(KIND=JPRB) ::&
 & 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(KIND=JPRB) ::&
 & ZCLDLD(KLON,KLEV)  , ZCLDLU(KLON,KLEV)&
 & , ZCLDSW(KLON,KLEV)  , ZCLD0(KLON,KLEV)&
 & , ZDT0(KLON)        &
 & , ZEMIS(KLON)        , ZEMIW(KLON)&
 & , ZFIWP(KLON)        , ZFLWP(KLON)      , ZFRWP(KLON)&
 & , ZIWC(KLON)         , ZLWC(KLON)&
 !cc            , ZRWC(KLON)
 & , ZMU0(KLON)         , ZOZ(KLON,KLEV)   , ZOZN(KLON,KLEV)&
 & , ZPMB(KLON,KLEV+1)  , ZPSOL(KLON)&
 & , ZTAVE (KLON,KLEV)  , ZTL(KLON,KLEV+1)&
 & , ZVIEW(KLON)  
REAL(KIND=JPRB) ::&
 & 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)&  
 & , ZDIRFS(KLON,NSW)   , ZDIFFS(KLON,NSW)  
REAL(KIND=JPRB) ::&
 & ZALFICE(KLON)      , ZGAMICE(KLON)     , ZBICE(KLON)   , ZDESR(KLON)&
 & , ZRADIP(KLON)       , ZRADLP(KLON)     &
 !cc           , ZRADRD(KLON)
 & , ZRAINT(KLON)       , ZRES(KLON)&
 & , ZTICE(KLON)        , ZEMIT(KLON),  ZBICFU(KLON)&
 & , ZKICFU(KLON)&
 & , PTOAG(KLON,JPGPT),PTOACG(KLON,JPGPT) !FC


REAL(KIND=JPRB) :: ZSUDU(KLON)   , ZPARF(KLON)       , ZUVDF(KLON), ZPARCF(KLON)
INTEGER(KIND=JPIM) :: IKL, JK, JKL, JKLP1, JKP1, JL, JNU, JRTM, JSW, INDLAY

INTEGER(KIND=JPIM) :: JLW , JI , J1, J2, J  !FC


REAL(KIND=JPRB) :: ZASYMX, ZDIFFD, ZGI, ZGL, ZGR, ZIWGKG, ZLWGKG,&
 & ZMSAID, ZMSAIU, ZMSALD, ZMSALU, ZRSAIA, ZRSAID, ZRSAIE, ZRSAIF, ZRSAIG, ZRSALD, &
 & ZMULTI, ZMULTL, ZOI   , ZOL, &
 & ZOMGMX, ZOR, ZRMUZ, ZRWGKG, ZTAUD, ZTAUMX, ZTEMPC, &
 & ZTOI, ZTOL, ZTOR, ZZFIWP, ZZFLWP, ZDPOG, ZPODT  

REAL(KIND=JPRB) :: ZALND, ZASEA, ZD, ZDEN, ZNTOT, ZNUM, ZRATIO, Z1RADI, &
 & Z1RADL, ZBETAI, ZOMGI, ZOMGP, ZFDEL, ZTCELS, ZFSR, ZAIWC, &
 & ZBIWC, ZTBLAY, ZADDPLK, ZPLANCK, ZEXTCF, Z1MOMG, &
 & ZDefRe, ZRefDe, ZVI , ZMABSD 

!REAL(KIND=JPRB) :: ZAVDP(KLON), ZAVTO(KLON), ZSQTO(KLON)
REAL(KIND=JPRB) :: ZAVTO(KLON), ZSQTO(KLON)
REAL(KIND=JPRB) :: ZSQUAR(KLON,KLEV), ZVARIA(KLON,KLEV)
INTEGER(KIND=JPIM) :: IKI, JKI, JEXPLR, JXPLDN
LOGICAL         :: LLDEBUG


REAL(KIND=JPRB) :: ZHOOK_HANDLE

#include "lw.intfb.h"
#include "rrtm_rrtm_140gp.intfb.h"
#include "sw.intfb.h"

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

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

IF (LHOOK) CALL DR_HOOK('RADLSW',0,ZHOOK_HANDLE)

LLDEBUG=.FALSE.
ZRefDe = RRe2De
ZDefRe = 1.0_JPRB / ZRefDe

DO JL = KIDIA,KFDIA
  ZFCUP(JL,KLEV+1) = 0.0_JPRB
  ZFCDWN(JL,KLEV+1) = REPLOG
  ZFSUP(JL,KLEV+1) = 0.0_JPRB
  ZFSDWN(JL,KLEV+1) = REPLOG
  PFLUX(JL,1,KLEV+1) = 0.0_JPRB
  PFLUX(JL,2,KLEV+1) = 0.0_JPRB
  PFLUC(JL,1,KLEV+1) = 0.0_JPRB
  PFLUC(JL,2,KLEV+1) = 0.0_JPRB
!FC
  DO JI = 1, JPGPT
  PTOAG(JL,JI) = 0.0_JPRB
  PTOACG(JL,JI) = 0.0_JPRB
  ENDDO
  DO JI = 1, NLW
  PTOAB(JL,JI) = 0.0_JPRB
  PTOACB(JL,JI) = 0.0_JPRB
  ENDDO
!FC

  ZFSDNN(JL) = 0.0_JPRB
  ZFSDNV(JL) = 0.0_JPRB
  ZFCDNN(JL) = 0.0_JPRB
  ZFCDNV(JL) = 0.0_JPRB
  ZFSUPN(JL) = 0.0_JPRB
  ZFSUPV(JL) = 0.0_JPRB
  ZFCUPN(JL) = 0.0_JPRB
  ZFCUPV(JL) = 0.0_JPRB
  ZPSOL(JL) = PAPH(JL,KLEV+1)
  ZPMB(JL,1) = ZPSOL(JL) / 100.0_JPRB
  ZDT0(JL) = PTS(JL) - PTH(JL,KLEV+1)
  PSUDU(JL) = 0.0_JPRB
  PPARF(JL) = 0.0_JPRB
  PPARCF(JL)= 0.0_JPRB
  PUVDF(JL) = 0.0_JPRB
  PSFSWDIR(JL,:)=0.0_JPRB
  PSFSWDIF(JL,:)=0.0_JPRB
  IBAS(JL) = INT ( 0.01_JPRB + PNBAS(JL) )
  ITOP(JL) = INT ( 0.01_JPRB + PNTOP(JL) )
ENDDO

!*         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)
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.0_JPRB

!-- ZOZ in cm.atm for SW scheme    
    ZOZ(JL,JK)   = POZON(JL,JKL) * 46.6968_JPRB / RG

    ZCLD0(JL,JK) = 0.0_JPRB
    ZFCUP(JL,JK) = 0.0_JPRB
    ZFCDWN(JL,JK) = 0.0_JPRB
    ZFSUP(JL,JK) = 0.0_JPRB
    ZFSDWN(JL,JK) = 0.0_JPRB
    PFLUX(JL,1,JK) = 0.0_JPRB
    PFLUX(JL,2,JK) = 0.0_JPRB
    PFLUC(JL,1,JK) = 0.0_JPRB
    PFLUC(JL,2,JK) = 0.0_JPRB
  ENDDO
ENDDO

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.0_JPRB
ENDDO
!***

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

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

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)  = 0.0_JPRB
      ZOMEGA(JL,JSW,JK)= 1.0_JPRB
      ZCG(JL,JSW,JK)   = 0.0_JPRB
    ENDDO
  ENDDO
  DO JL = KIDIA,KFDIA
    ZCLDSW(JL,JK)  = 0.0_JPRB
    ZCLDLD(JL,JK)  = 0.0_JPRB
    ZCLDLU(JL,JK)  = 0.0_JPRB
  ENDDO

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

  DO JL = KIDIA,KFDIA

! --- LIQUID WATER CONTENT (g.m-3) AND LIQUID WATER PATH (g.m-2)
    IF (PCLFR(JL,IKL) > REPSC ) THEN
      ZLWGKG=MAX(PQLWP(JL,IKL)*1000.0_JPRB,0.0_JPRB)
      ZIWGKG=MAX(PQIWP(JL,IKL)*1000.0_JPRB,0.0_JPRB)
      ZLWGKG=ZLWGKG/PCLFR(JL,IKL)
      ZIWGKG=ZIWGKG/PCLFR(JL,IKL)
    ELSE
      ZLWGKG=0.0_JPRB
      ZIWGKG=0.0_JPRB
    ENDIF
    ZRWGKG=0.0_JPRB
    ZRAINT(JL)=0.0_JPRB

! --- RAIN LIQUID WATER CONTENT (g.m-3) AND LIQUID WATER PATH (g.m-2)
!    IF (PRAINT(JL,IKL) >= REPSCW) 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
!    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

  ENDDO
  DO JL = KIDIA,KFDIA
! --- EFFECTIVE RADIUS FOR WATER, ICE AND RAIN PARTICLES

! very old parametrization as f(pressure)

    IF (NRADLP == 0) THEN
!-- very old parametrization as f(pressure) ERA-15
      ZRADLP(JL)=10.0_JPRB + (100000.0_JPRB-PAP(JL,IKL))*3.5_JPRB

    ELSEIF (NRADLP == 1) THEN
! simple distinction between land (10) and ocean (13) Zhang and Rossow
      IF (PLSM(JL) < 0.5_JPRB) THEN
        ZRADLP(JL)=13.0_JPRB
      ELSE
        ZRADLP(JL)=10.0_JPRB
      ENDIF
      
    ELSEIF (NRADLP == 2) THEN
!--  based on Martin et al., 1994, JAS
      IF (PLSM(JL) < 0.5_JPRB) THEN
        IF (LCCNO) THEN
!          ZASEA=50.0_JPRB
          ZASEA=PCCNO(JL)
        ELSE  
          ZASEA=RCCNSEA
        ENDIF  
        ZD=0.33_JPRB
        ZNTOT=-1.15E-03_JPRB*ZASEA*ZASEA+0.963_JPRB*ZASEA+5.30_JPRB
      ELSE
        IF (LCCNL) THEN 
!          ZALND=900.0_JPRB
          ZALND=PCCNL(JL)
        ELSE  
          ZALND=RCCNLND
        ENDIF  
        ZD=0.43_JPRB
        ZNTOT=-2.10E-04_JPRB*ZALND*ZALND+0.568_JPRB*ZALND-27.9_JPRB
      ENDIF
      ZNUM=3.0_JPRB*ZLWC(JL)*(1.0_JPRB+3.0_JPRB*ZD*ZD)**2
      ZDEN=4.0_JPRB*RPI*ZNTOT*(1.0_JPRB+ZD*ZD)**3
      IF((ZNUM/ZDEN) > REPLOG)THEN
        ZRADLP(JL)=100.0_JPRB*EXP(0.333_JPRB*LOG(ZNUM/ZDEN))
        ZRADLP(JL)=MAX(ZRADLP(JL), 4.0_JPRB)
        ZRADLP(JL)=MIN(ZRADLP(JL),16.0_JPRB)
      ELSE
        ZRADLP(JL)=4.0_JPRB
      ENDIF

    ELSEIF (NRADLP == 3) THEN  
! one uses the cloud droplet radius from newmicro
! IKL or JK ?? - I think IKL but needs to be verified > ref_liq_i
! (inverted) is used in the call of RECMWF_AERO in radlwsw_m.F90,
! so everything is fine - JBM 6/2019
        ZRADLP(JL)=PREF_LIQ(JL,IKL)
    ENDIF  

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

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

  ENDDO
  DO JL = KIDIA,KFDIA

! diagnosing the ice particle effective radius/diameter

!- 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))    

    IF (NRADIP == 0) THEN
!-- fixed 40 micron effective radius
      ZRADIP(JL)= 40.0_JPRB
      ZDESR(JL) = ZDefRe * ZRADIP(JL)
      
    ELSEIF (NRADIP == 1) THEN 

!-- old formulation based on Liou & Ou (1994) temperature (40-130microns)    
      ZRADIP(JL)=MAX(ZRADIP(JL),40.0_JPRB)
      ZDESR(JL) = ZDefRe * ZRADIP(JL)
      
    ELSEIF (NRADIP == 2) THEN  
!-- formulation following Jakob, Klein modifications to ice content    
      ZRADIP(JL)=MAX(ZRADIP(JL),30.0_JPRB)
      ZRADIP(JL)=MIN(ZRADIP(JL),60.0_JPRB)
      ZDESR(JL)= ZDefRe * ZRADIP(JL)
 
    ELSEIF (NRADIP == 3  ) THEN
 
!- ice particle effective radius =f(T,IWC) from Sun and Rikus (1999)
! revised by Sun (2001)
      IF (ZIWC(JL) > 0.0_JPRB ) THEN
        ZTEMPC = PT(JL,IKL)-83.15_JPRB
        ZTCELS = PT(JL,IKL)-RTT
        ZFSR = 1.2351_JPRB +0.0105_JPRB * ZTCELS
! Sun, 2001 (corrected from Sun & Rikus, 1999)
        ZAIWC = 45.8966_JPRB * ZIWC(JL)**0.2214_JPRB
        ZBIWC = 0.7957_JPRB * ZIWC(JL)**0.2535_JPRB
        ZDESR(JL) = ZFSR * (ZAIWC + ZBIWC*ZTEMPC)
!-new        ZDESR(JL) = MIN ( MAX( ZDESR(JL), 30.0_JPRB), 155.0_JPRB)
        ZDESR(JL) = MIN ( MAX( ZDESR(JL), 45.0_JPRB), 350.0_JPRB)
        ZRADIP(JL)= ZRefDe * ZDESR(JL)
      ELSE
!        ZDESR(JL) = 92.5_JPRB
        ZDESR(JL) = 80.0_JPRB
        ZRADIP(JL)= ZRefDe * ZDESR(JL)
      ENDIF  
 
    ELSEIF (NRADIP == 4  ) THEN
! one uses the cloud droplet radius from newmicro
! IKL or JK ?? - I think IKL but needs to be verified
        ZRADIP(JL)=PREF_ICE(JL,IKL)
    ENDIF  
    
  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=0.0_JPRB
      ZGL =0.0_JPRB
      ZOL =0.0_JPRB
      ZTOI=0.0_JPRB
      ZGI =0.0_JPRB
      ZOI =0.0_JPRB
      ZTOR=0.0_JPRB
      ZGR =0.0_JPRB
      ZOR =0.0_JPRB
      IF (ZFLWP(JL)+ZFIWP(JL)+ZFRWP(JL) > 2.0_JPRB * REPSCW ) THEN
        IF (ZFLWP(JL) >= REPSCW ) THEN
          IF (NLIQOPT /= 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)
!-- NB: RSWINHF is there simply for making the CY29R2 branch bit compatible with 
! the previous. Should be cleaned when RRTM_SW becomes active
            ZOL  = RYFWCC(JSW)-RYFWCD(JSW)*EXP(-RYFWCE(JSW)*ZTOL*RSWINHF)
          ENDIF 
        ENDIF

        IF (ZFIWP(JL) >= REPSCW ) THEN
          IF (NICEOPT <= 1) THEN
!-- SW: Ebert-Curry          
            ZTOI = ZFIWP(JL)*(REBCUA(JSW)+REBCUB(JSW)/ZRADIP(JL))
            ZGI  = REBCUE(JSW)+REBCUF(JSW)*ZRADIP(JL)
            ZOI  = 1.0_JPRB - REBCUC(JSW)-REBCUD(JSW)*ZRADIP(JL)
            
          ELSEIF (NICEOPT == 2) THEN  
!-- SW: Fu-Liou 1993
            Z1RADI = 1.0_JPRB / ZDESR(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  = 1.0_JPRB - 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.0_JPRB -ZFDEL)*ZOMGP + ZFDEL*3.0_JPRB) / 3.0_JPRB
            
          ELSEIF (NICEOPT == 3) THEN  
!-- SW: Fu 1996
            Z1RADI = 1.0_JPRB / 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  = 1.0_JPRB - ZOMGI
            ZGI  = RFUCC0(JSW)+ZDESR(JL)*(RFUCC1(JSW) + ZDESR(JL) &
             &   *(RFUCC2(JSW)+ZDESR(JL)* RFUCC3(JSW) )) 
            ZGI  = MIN(1.0_JPRB, ZGI)
      
          ENDIF
        ENDIF

!        IF (ZFRWP(JL) >= REPSCW ) THEN
!          ZTOR= ZFRWP(JL)*0.003_JPRB * ZRAINT(JL)**(-0.22_JPRB)         
!          ZOR = 1.0_JPRB - RROMA(JSW)*ZRAINT(JL)**RROMB(JSW)
!          ZGR = RRASY(JSW)
!        ENDIF   

!  - 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

! --- 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
    ENDDO
  ENDDO

  IF(LLDEBUG) THEN
   call writefield_phy("radlsw_ztau",ztau(:,1,:),klev)
  ENDIF

!          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)=0.0_JPRB
      ZGAMICE(JL)=0.0_JPRB
      ZBICE(JL)=0.0_JPRB
      ZTICE(JL)=(PT(JL,IKL)-TSTAND)/TSTAND
      IF (NICEOPT == 1) THEN
        ZBICFU(JL)=1.0_JPRB
      ELSE
        ZBICFU(JL)=0.0_JPRB
      ENDIF
      ZKICFU(JL)=0.0_JPRB
    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    
    IF (NICEOPT == 2 .OR. NICEOPT == 3) THEN
      DO JRTM=1,16
        DO JL=KIDIA,KFDIA
          IF (PT(JL,IKL) < 160.0_JPRB) THEN
            INDLAY=1
            ZTBLAY =PT(JL,IKL)-160.0_JPRB
          ELSEIF (PT(JL,IKL) < 339.0_JPRB ) THEN
            INDLAY=PT(JL,IKL)-159.0_JPRB
            INDLAY=MAX(INDLAY,1)
            ZTBLAY =PT(JL,IKL)-INT(PT(JL,IKL))
          ELSE 
            INDLAY=180
            ZTBLAY =PT(JL,IKL)-339.0_JPRB
          ENDIF
          ZADDPLK = TOTPLNK(INDLAY+1,JRTM)-TOTPLNK(INDLAY,JRTM)
          ZPLANCK = DELWAVE(JRTM) * (TOTPLNK(INDLAY,JRTM) + ZTBLAY*ZADDPLK)
          ZBICFU(JL) = ZBICFU(JL) + ZPLANCK
        
          IF (ZIWC(JL) > 0.0_JPRB ) THEN
            ZRATIO =  1.0_JPRB / ZDESR(JL) 
            IF (NICEOPT == 2) THEN
! ice cloud spectral emissivity a la Fu & Liou (1993)
              ZMABSD = RFULIO(JRTM,1) + ZRATIO &
               & *(RFULIO(JRTM,2) + ZRATIO*RFULIO(JRTM,3))  
          
! ice cloud spectral emissivity a la Fu et al (1998)
            ELSEIF (NICEOPT == 3) THEN 
              ZMABSD = RFUETA(JRTM,1) + ZRATIO &
               & *(RFUETA(JRTM,2) + ZRATIO*RFUETA(JRTM,3))  
            ENDIF
            ZKICFU(JL) = ZKICFU(JL)+ ZMABSD*ZPLANCK
          ENDIF  
        ENDDO
      ENDDO
    ENDIF
    
    DO JL = KIDIA,KFDIA
      ZGAMICE(JL) = ZGAMICE(JL) / ZBICE(JL)
      ZALFICE(JL) = ZALFICE(JL) / ZBICE(JL)
      ZKICFU(JL)  = ZKICFU(JL) / ZBICFU(JL)
      
      IF (ZFLWP(JL)+ZFIWP(JL) > REPSCW) THEN

        IF (NLIQOPT == 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
          
        ENDIF  
          
        IF (NICEOPT == 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

        ELSEIF (NICEOPT == 1) THEN
! ice cloud emissivity a la Ebert & Curry (1992)
          ZMSAID= 1.66_JPRB*(ZALFICE(JL)+ZGAMICE(JL)/ZRADIP(JL))
          ZMSAIU= ZMSAID
          
        ELSEIF (NICEOPT == 2 .OR. NICEOPT == 3) THEN  
! ice cloud emissivity a la Fu & Liou (1993) or Fu et al. (1998)
          ZMSAID= 1.66_JPRB*ZKICFU(JL)
          ZMSAIU= ZMSAID          
        ENDIF 
       
        IF (NINHOM == 1) THEN
          ZZFLWP= ZFLWP(JL) * RLWINHF
          ZZFIWP= ZFIWP(JL) * RLWINHF
        ELSE
          ZZFLWP= ZFLWP(JL)
          ZZFIWP= ZFIWP(JL)
        ENDIF

! effective cloudiness accounting for condensed water
        ZCLDLD(JL,JK) = PCLFR(JL,IKL)*(1.0_JPRB-EXP(-ZMSALD*ZZFLWP-ZMSAID* &
         & ZZFIWP))  
        ZCLDLU(JL,JK) = PCLFR(JL,IKL)*(1.0_JPRB-EXP(-ZMSALU*ZZFLWP-ZMSAIU* &
         & ZZFIWP))  
      ENDIF
    ENDDO

  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

    DO JRTM=1,16
      DO JL = KIDIA,KFDIA
        ZTAUCLD(JL,JK,JRTM) = 0.0_JPRB
        ZMSALD = 0.0_JPRB
        ZMSAID = 0.0_JPRB
        
        IF (ZFLWP(JL)+ZFIWP(JL) > REPSCW) THEN
    
          IF (NLIQOPT == 0 .OR. NLIQOPT >= 3 ) THEN
! water cloud total emissivity a la Smith and Shi (1992)
            ZMULTL=1.2_JPRB-0.006_JPRB*ZRADLP(JL)
            ZRSALD= 0.144_JPRB*ZMULTL / 1.66_JPRB
            
          ELSEIF (NLIQOPT == 1) THEN
! water cloud spectral emissivity a la Savijarvi (1997)
            ZRSALD= RHSAVI(JRTM,1) + ZRADLP(JL)&
             & *(RHSAVI(JRTM,2) + ZRADLP(JL)*RHSAVI(JRTM,3))  
             
          ELSEIF (NLIQOPT == 2) THEN
! water cloud spectral emissivity a la Lindner and Li (2000)
            Z1RADL = 1.0_JPRB / ZRADLP(JL)
            ZEXTCF = RLILIA(JRTM,1)+ZRADLP(JL)*RLILIA(JRTM,2)+ Z1RADL*&
             & (RLILIA(JRTM,3) + Z1RADL*(RLILIA(JRTM,4) + Z1RADL*&
             & RLILIA(JRTM,5) ))  
            Z1MOMG = RLILIB(JRTM,1) + Z1RADL*RLILIB(JRTM,2) &
             & + ZRADLP(JL) *(RLILIB(JRTM,3) + ZRADLP(JL)*RLILIB(JRTM,4) )
            ZRSALD = Z1MOMG * ZEXTCF
          ENDIF  
         
          IF (NICEOPT == 0) THEN
! ice cloud spectral emissivity a la Smith & Shi (1992)
            ZMULTI=1.2_JPRB-0.006_JPRB*ZRADIP(JL)
            ZRSAID= 0.103_JPRB*ZMULTI / 1.66_JPRB
            
          ELSEIF (NICEOPT == 1) THEN
! ice cloud spectral emissivity a la Ebert-Curry (1992)
            ZRSAID= REBCUH(JRTM)+REBCUG(JRTM)/ZRADIP(JL)
            
          ELSEIF (NICEOPT == 2) THEN
! ice cloud spectral emissivity a la Fu & Liou (1993)
            Z1RADI = 1.0_JPRB / ZDESR(JL)
            ZRSAID = RFULIO(JRTM,1) + Z1RADI &
             & *(RFULIO(JRTM,2) + Z1RADI * RFULIO(JRTM,3))  
             
          ELSEIF (NICEOPT == 3) THEN
! ice cloud spectral emissivity a la Fu et al (1998) including 
! parametrisation for LW scattering effect  
            Z1RADI = 1.0_JPRB / ZDESR(JL)
            ZRSAIE = RFUETA(JRTM,1) + Z1RADI &
             &*(RFUETA(JRTM,2) + Z1RADI * RFUETA(JRTM,3)) 
            ZRSAIA = Z1RADI*(RFUETB(JRTM,1) +ZDESR(JL)*( RFUETB(JRTM,2) +ZDESR(JL)*( RFUETB(JRTM,3) +ZDESR(JL)* RFUETB(JRTM,4))))
            ZRSAIG = RFUETC(JRTM,1) +ZDESR(JL)*( RFUETC(JRTM,2) +ZDESR(JL)*( RFUETC(JRTM,3) +ZDESR(JL)* RFUETC(JRTM,4))) 
            ZRSAIF = 0.5_JPRB + ZRSAIG*( 0.3738_JPRB + ZRSAIG*( 0.0076_JPRB + ZRSAIG*0.1186_JPRB ) )
            ZRSAID = (1.0_JPRB - ZRSAIA/ZRSAIE * ZRSAIF) * ZRSAIE
          ENDIF    
         
          ZTAUD = ZRSALD*ZFLWP(JL)+ZRSAID*ZFIWP(JL)

! Diffusivity correction within clouds a la Savijarvi
          IF (LDIFFC) THEN
            ZDIFFD=MIN(MAX(1.517_JPRB-0.156_JPRB*LOG(ZTAUD) , 1.0_JPRB), &
             &     2.0_JPRB)
          ELSE
            ZDIFFD=1.66_JPRB
          ENDIF

          ZTAUCLD(JL,JK,JRTM) = ZTAUD*ZDIFFD
        ENDIF
        
      ENDDO
    ENDDO

  ENDIF

ENDDO

NUAER = NUA
NTRAER = NTRA

!     ------------------------------------------------------------------
! 
!          2.6    SCALING OF OPTICAL THICKNESS
!                 SPECTRALLY, ACCOUNTING FOR VERTICAL VARIABILITY

JEXPLR=NLAYINH
JXPLDN=2*JEXPLR+1

IF (NINHOM == 1) THEN
!-- simple scaling a la Tiedtke (1996) with RSWINHF in SW and RLWINHF in LW
  DO JSW=1,NSW
    DO JK=1,KLEV
      DO JL=KIDIA,KFDIA
        ZTAU(JL,JSW,JK)=ZTAU(JL,JSW,JK) * RSWINHF
      ENDDO
    ENDDO
  ENDDO

  DO JRTM=1,16
    DO JK=1,KLEV
      DO JL=KIDIA,KFDIA
        ZTAUCLD(JL,JK,JRTM)=ZTAUCLD(JL,JK,JRTM) * RLWINHF
      ENDDO
    ENDDO
  ENDDO

ELSEIF (JEXPLR /= 0) THEN
  DO JSW=1,NSW
    DO JK=1,KLEV
      DO JL=KIDIA,KFDIA
        ZSQUAR(JL,JK)=0.0_JPRB
        ZVARIA(JL,JK)=1.0_JPRB
      ENDDO
    ENDDO
!-- range should be defined from Hogan & Illingworth
    DO JK=1+JEXPLR,KLEV-JEXPLR
      DO JL=KIDIA,KFDIA
!        ZAVDP(JL)=0.0_JPRB
        ZAVTO(JL)=0.0_JPRB
        ZSQTO(JL)=0.0_JPRB
      ENDDO
      DO JKI=JK-JEXPLR,JK+JEXPLR
        IKI=KLEV+1-JKI
        DO JL=KIDIA,KFDIA
!          ZAVDP(JL)=ZAVDP(JL)+PDP(JL,IKI)/RG
          ZAVTO(JL)=ZAVTO(JL)+ZTAU(JL,JSW,JKI)
        ENDDO
      ENDDO
      DO JL=KIDIA,KFDIA
!        ZAVTO(JL)=ZAVTO(JL)/ZAVDP(JL)
        ZAVTO(JL)=ZAVTO(JL)/JXPLDN
      ENDDO
      DO JKI=JK-JEXPLR,JK+JEXPLR
        IKI=KLEV+1-JKI
        DO JL=KIDIA,KFDIA
!          ZSQTO(JL)=ZSQTO(JL)+(ZTAU(JL,JSW,JKI)/PDP(JL,IKI)-ZAVTO(JL))**2
          ZSQTO(JL)=ZSQTO(JL)+(ZTAU(JL,JSW,JKI)-ZAVTO(JL))**2
        ENDDO
      ENDDO
      DO JL=KIDIA,KFDIA
        ZSQTO(JL)=SQRT(ZSQTO(JL)/(JXPLDN*(JXPLDN-1)))
        IF (ZAVTO(JL) > 0.0_JPRB) THEN
          ZVARIA(JL,JK)=(ZSQTO(JL)/ZAVTO(JL))**2
          ZSQUAR(JL,JK)=EXP(-ZVARIA(JL,JK))
        ELSE
          ZVARIA(JL,JK)=0.0_JPRB
          ZSQUAR(JL,JK)=1.0_JPRB
        ENDIF

!-- scaling a la Barker
        IF (NINHOM ==2) THEN
          ZTAU(JL,JSW,JK)=ZTAU(JL,JSW,JK)*ZSQUAR(JL,JK)

!-- scaling a la Cairns et al.
        ELSEIF (NINHOM == 3) THEN
          ZVI=ZVARIA(JL,JK) 
          ZTAU(JL,JSW,JK)  = ZTAU(JL,JSW,JK)/(1.0_JPRB+ZVI)
          ZOMEGA(JL,JSW,JK)= ZOMEGA(JL,JSW,JK) &
            &   /(1.0_JPRB + ZVI*(1.0_JPRB-ZOMEGA(JL,JSW,JK) ) )
          ZCG(JL,JSW,JK)   = ZCG(JL,JSW,JK) &
            & *(1.0_JPRB+ZVI*(1.0_JPRB-ZOMEGA(JL,JSW,JK))) &
            & /(1.0_JPRB+ZVI*(1.0_JPRB-ZOMEGA(JL,JSW,JK)*ZCG(JL,JSW,JK)))
        ENDIF
      ENDDO
!      JL=KIDIA
!      print 9261,JSW,JK,ZTAU(JL,JSW,JK),ZAVTO(JL),ZSQTO(JL),ZVARIA(JL,JK),ZSQUAR(JL,JK)
9261   format(1x,'Varia1 ',2I3,7F10.4)
    ENDDO
  ENDDO


  DO JRTM=1,16
    DO JK=1,KLEV
      DO JL=KIDIA,KFDIA
        ZSQUAR(JL,JK)=0.0_JPRB
        ZVARIA(JL,JK)=1.0_JPRB
      ENDDO
    ENDDO
!-- range to be defined from Hogan & Illingworth
    DO JK=1+JEXPLR,KLEV-JEXPLR
      DO JL=KIDIA,KFDIA
!        ZAVDP(JL)=0.0_JPRB
        ZAVTO(JL)=0.0_JPRB
        ZSQTO(JL)=0.0_JPRB
      ENDDO
      DO JKI=JK-JEXPLR,JK+JEXPLR
        IKI=KLEV+1-JKI
        DO JL=KIDIA,KFDIA
!          ZAVDP(JL)=ZAVDP(JL)+PDP(JL,IKI)/RG
          ZAVTO(JL)=ZAVTO(JL)+ZTAUCLD(JL,JKI,JRTM)
        ENDDO
      ENDDO
      DO JL=KIDIA,KFDIA
!        ZAVTO(JL)=ZAVTO(JL)/ZAVDP(JL)
        ZAVTO(JL)=ZAVTO(JL)/JXPLDN
      ENDDO
      DO JKI=JK-JEXPLR,JK+JEXPLR
        IKI=KLEV+1-JKI
        DO JL=KIDIA,KFDIA
!          ZSQTO(JL)=ZSQTO(JL)+(ZTAUCLD(JL,JKI,JRTM)/PDP(JL,IKI)-ZAVTO(JL))**2
            ZSQTO(JL)=ZSQTO(JL)+(ZTAUCLD(JL,JKI,JRTM)-ZAVTO(JL))**2
        ENDDO
      ENDDO
      DO JL=KIDIA,KFDIA
        ZSQTO(JL)=SQRT(ZSQTO(JL)/(JXPLDN*(JXPLDN-1)))
        IF (ZAVTO(JL) > 0.0_JPRB) THEN
          ZVARIA(JL,JK)=(ZSQTO(JL)/ZAVTO(JL))**2
          ZSQUAR(JL,JK)=EXP(-ZVARIA(JL,JK))
        ELSE
          ZVARIA(JL,JK)=0.0_JPRB
          ZSQUAR(JL,JK)=1.0_JPRB
        ENDIF

!-- scaling a la Barker
        IF (NINHOM ==2) THEN
          ZTAUCLD(JL,JK,JRTM)=ZTAUCLD(JL,JK,JRTM)*ZSQUAR(JL,JK)

!-- scaling a la Cairns et al.
        ELSEIF (NINHOM == 3) THEN
          ZVI=ZVARIA(JL,JK) 
          ZTAUCLD(JL,JK,JRTM)=ZTAUCLD(JL,JK,JRTM)/(1.0_JPRB+ZVI)
        ENDIF
      ENDDO
!      JL=KIDIA
!      print 9262,JRTM,JK,ZTAUCLD(JL,JK,JRTM),ZAVTO(JL),ZSQTO(JL),ZVARIA(JL,JK),ZSQUAR(JL,JK)
9262   format(1x,'Varia2 ',2I3,7F10.4)
    ENDDO
  ENDDO
ENDIF



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

DO JL = KIDIA,KFDIA
  ZVIEW(JL) = DIFF
ENDDO

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

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

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

!print *,'RADLSW: LPHYLIN, LRRTM',LPHYLIN, LRRTM
IF (.NOT.LPHYLIN) THEN
  IF ( .NOT. LRRTM) THEN
           !FC ON EST DANS LE CAS OU NOT.RRTM EST TRUE DONC ON NE PASSE PAS LA


    CALL LW &
     & ( KIDIA , KFDIA , KLON  , KLEV , KMODE,&
     & PCCO2 , ZCLDLD, ZCLDLU,&
     & PDP   , ZDT0  , ZEMIS , ZEMIW,&
     & ZPMB  , POZON , ZTL,&
     & PAER  , ZTAVE , ZVIEW , PQ,&
     & ZEMIT , PFLUX , PFLUC &
     & )  
!   print *,'RADLSW: apres CALL LW'
    IF(LLDEBUG) THEN
    call writefield_phy('radlsw_flux1',PFLUX(:,1,:),klev+1)
    call writefield_phy('radlsw_flux2',PFLUX(:,2,:),klev+1)
    call writefield_phy('radlsw_fluc1',PFLUC(:,1,:),klev+1)
    call writefield_phy('radlsw_fluc2',PFLUC(:,2,:),klev+1)
    ENDIF

  ELSE

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

!  i)  pass ZOZN (ozone mass mixing ratio) 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

!   print *,'RADLSW: avant CALL RRTM_RRTM_140GP,PAP=',PAP(1,:)
    CALL RRTM_RRTM_140GP &
     & ( KIDIA , KFDIA , KLON  , KLEV,&
     & PAER  , PAPH  , PAP,&
     & PTS   , PTH   , PT,&
     & ZEMIS , ZEMIW,&
     & PQ    , PCCO2 , ZOZN  ,&
     & ZCLDSW  , ZTAUCLD,&
     & PTAU_LW,&
     & ZEMIT , PFLUX , PFLUC , ZTCLEAR,&
     & PTOAG, PTOACG )  !FC 
!   print *,'RADLSW: apres CALL RRTM_RRTM_140GP'

  ENDIF
ELSE
  ZEMIT (:)   = 0.0_JPRB
  PFLUX(:,:,:)= 0.0_JPRB
  PFLUC(:,:,:)= 0.0_JPRB
! print *,'RADLSW: ZEMIT,PFLUX et PFLUC = 0'
ENDIF

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

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

ZRMUZ=0.0_JPRB
DO JL = KIDIA,KFDIA
  ZRMUZ = MAX (ZRMUZ, ZMU0(JL))
ENDDO

IF (NSTEP == 0 .AND. LEDBUG .AND. ZMU0(KIDIA) > 0.0_JPRB) THEN
  WRITE(NULOUT,'(4E15.8)') PRII0,PCCO2,ZPSOL(KIDIA),ZMU0(KIDIA)
  WRITE(NULOUT,'("ZALBD ",6E15.8)') (ZALBD(KIDIA,JSW),JSW=1,NSW)
  WRITE(NULOUT,'("ZALBP ",6E15.8)') (ZALBP(KIDIA,JSW),JSW=1,NSW)
  WRITE(NULOUT,'("PQ    ",10E12.5)') (PQ(KIDIA,JK),JK=1,KLEV)
  WRITE(NULOUT,'("PQS   ",10E12.5)') (PQS(KIDIA,JK),JK=1,KLEV)
  WRITE(NULOUT,'("PDP   ",10E12.5)') (PDP(KIDIA,JK),JK=1,KLEV)
  WRITE(NULOUT,'("ZPMB  ",10E12.5)') (ZPMB(KIDIA,JK),JK=1,KLEV+1)
  WRITE(NULOUT,'("ZTAVE ",10E12.5)') (ZTAVE(KIDIA,JK),JK=1,KLEV)
  WRITE(NULOUT,'("ZCLDSW",10E12.5)') (ZCLDSW(KIDIA,JK),JK=1,KLEV)
  WRITE(NULOUT,'("ZTAU  ",10E12.5)') ((ZTAU(KIDIA,JSW,JK),JK=1,KLEV),JSW=1,NSW)
  WRITE(NULOUT,'("ZCG   ",10E12.5)') ((ZCG(KIDIA,JSW,JK),JK=1,KLEV),JSW=1,NSW)
  WRITE(NULOUT,'("ZOMEGA",10E12.5)') ((ZOMEGA(KIDIA,JSW,JK),JK=1,KLEV),JSW=1,NSW)
  WRITE(NULOUT,'("ZOZ   ",10E12.5)') (ZOZ(KIDIA,JK),JK=1,KLEV)
  WRITE(NULOUT,'("PAER  ",10E12.5)') ((PAER(KIDIA,JSW,JK),JK=1,KLEV),JSW=1,NSW)
ENDIF

IF (NSTEP == 0 .AND. LEDBUG .AND. ZMU0(KIDIA) > 0.0_JPRB) THEN
  WRITE(NULOUT,'(4E15.8)') PRII0,PCCO2,ZPSOL(KIDIA),ZMU0(KIDIA)
  WRITE(NULOUT,'("ZALBD ",6E15.8)') (ZALBD(KIDIA,JSW),JSW=1,NSW)
  WRITE(NULOUT,'("ZALBP ",6E15.8)') (ZALBP(KIDIA,JSW),JSW=1,NSW)
  WRITE(NULOUT,'("PQ    ",10E12.5)') (PQ(KIDIA,JK),JK=1,KLEV)
  WRITE(NULOUT,'("PQS   ",10E12.5)') (PQS(KIDIA,JK),JK=1,KLEV)
  WRITE(NULOUT,'("PDP   ",10E12.5)') (PDP(KIDIA,JK),JK=1,KLEV)
  WRITE(NULOUT,'("ZPMB  ",10E12.5)') (ZPMB(KIDIA,JK),JK=1,KLEV+1)
  WRITE(NULOUT,'("ZTAVE ",10E12.5)') (ZTAVE(KIDIA,JK),JK=1,KLEV)
  WRITE(NULOUT,'("ZCLDSW",10E12.5)') (ZCLDSW(KIDIA,JK),JK=1,KLEV)
  WRITE(NULOUT,'("ZTAU  ",10E12.5)') ((ZTAU(KIDIA,JSW,JK),JK=1,KLEV),JSW=1,NSW)
  WRITE(NULOUT,'("ZCG   ",10E12.5)') ((ZCG(KIDIA,JSW,JK),JK=1,KLEV),JSW=1,NSW)
  WRITE(NULOUT,'("ZOMEGA",10E12.5)') ((ZOMEGA(KIDIA,JSW,JK),JK=1,KLEV),JSW=1,NSW)
  WRITE(NULOUT,'("ZOZ   ",10E12.5)') (ZOZ(KIDIA,JK),JK=1,KLEV)
  WRITE(NULOUT,'("PAER  ",10E12.5)') ((PAER(KIDIA,JSW,JK),JK=1,KLEV),JSW=1,NSW)
ENDIF
CALL SW &
 & ( KIDIA , KFDIA , KLON  , KLEV  , KAER,&
 & PRII0 , PCCO2 , ZPSOL , ZALBD , ZALBP , PQ   , PQS,&
 & ZMU0  , ZCG   , ZCLDSW, PDP   , ZOMEGA, ZOZ  , ZPMB,&
 & ZTAU  , ZTAVE , PAER,&
 & PFSDN , PFSUP , PFSCDN, PFSCUP,&
 & ZFSDNN, ZFSDNV, ZFSUPN, ZFSUPV,&
 & ZFCDNN, ZFCDNV, ZFCUPN, ZFCUPV,&
 & ZSUDU , ZUVDF , ZPARF ,ZPARCF, ZDIFFS, ZDIRFS, &
 & LRDUST,PPIZA_DST,PCGA_DST,PTAUREL_DST&
   & )
PFSDNV=ZFSDNV
PFSDNN=ZFSDNN
IF (SIZE(PSFSWDIR,2)>1) THEN
  PSFSWDIR= ZDIRFS
  PSFSWDIF= ZDIFFS
ELSE
  PSFSWDIR (:,1) = ZFSDNV(:) + ZFSDNN(:)
  PSFSWDIF (:,:) = 0.
ENDIF

IF (NSTEP == 0 .AND. LEDBUG .AND. ZMU0(KIDIA) > 0.0_JPRB) THEN
  WRITE(NULOUT,'("ZFSDWN",10E12.5)') (ZFSDWN(KIDIA,JK),JK=1,KLEV)
  WRITE(NULOUT,'("ZFSUP ",10E12.5)') (ZFSUP (KIDIA,JK),JK=1,KLEV)
  WRITE(NULOUT,'("ZFCDWN",10E12.5)') (ZFCDWN(KIDIA,JK),JK=1,KLEV)
  WRITE(NULOUT,'("ZFCUP ",10E12.5)') (ZFCUP (KIDIA,JK),JK=1,KLEV)
  LEDBUG=.FALSE. 
ENDIF
IF (NSTEP == 0 .AND. LEDBUG .AND. ZMU0(KIDIA) > 0.0_JPRB) THEN
  WRITE(NULOUT,'("ZFSDWN",10E12.5)') (ZFSDWN(KIDIA,JK),JK=1,KLEV)
  WRITE(NULOUT,'("ZFSUP ",10E12.5)') (ZFSUP (KIDIA,JK),JK=1,KLEV)
  WRITE(NULOUT,'("ZFCDWN",10E12.5)') (ZFCDWN(KIDIA,JK),JK=1,KLEV)
  WRITE(NULOUT,'("ZFCUP ",10E12.5)') (ZFCUP (KIDIA,JK),JK=1,KLEV)
  LEDBUG=.FALSE. 
ENDIF
!     ------------------------------------------------------------------

!*         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
    PFLS(JL,JKL) = ZFSDWN(JL,JK) - ZFSUP(JL,JK)
    PFLT(JL,JKL) = - PFLUX(JL,1,JK) - PFLUX(JL,2,JK)
    PFCS(JL,JKL) = ZFCDWN(JL,JK) - ZFCUP(JL,JK)
    PFCT(JL,JKL) = - PFLUC(JL,1,JK) - PFLUC(JL,2,JK)
  ENDDO
ENDDO

!FC
  DO JL = KIDIA,KFDIA
    J1 = 1
    J2 = 0
    DO JI = 1, NLW  !NLW=16
      J2 = J2 + GP_PER_LWBAND(JI)
      DO J = J1, J2
        PTOAB(JL,JI) = PTOAB(JL,JI) + PTOAG(JL,J)
        PTOACB(JL,JI) = PTOACB(JL,JI) + PTOACG(JL,J)
      ENDDO
      J1 = J1 + GP_PER_LWBAND(JI)
    ENDDO
  ENDDO
!FC

DO JL = KIDIA,KFDIA
  PFRSOD(JL)=ZFSDWN(JL,1)
  PEMIT (JL)=ZEMIT (JL)
  PSUDU (JL)=ZSUDU (JL)
  PUVDF (JL)=ZUVDF (JL)
  PPARF (JL)=ZPARF (JL)
  PPARCF(JL)=ZPARCF(JL)
  PTINCF(JL)=PRII0 * ZMU0(JL) 
ENDDO
!print 9501,(PUVDF(JL),JL=KIDIA,KFDIA)
9501 format(1x,'RADLSW PUVDF: ',30f6.1)
!print 9502,(PPARF(JL),JL=KIDIA,KFDIA)
9502 format(1x,'RADLSW PPARF: ',30f6.1)

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

IF (LHOOK) CALL DR_HOOK('RADLSW',1,ZHOOK_HANDLE)
END SUBROUTINE RADLSW


































