radlsw.F90

Go to the documentation of this file.

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)

use write_field_phy

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

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)   :: 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)
REAL(KIND=JPRB) :: ZSUDU(klon)   , ZPARF(klon)       , ZUVDF(klon), ZPARCF(klon)
INTEGER(KIND=JPIM) :: IKL, JK, JKL, JKLP1, JKP1, JL, JNU, JRTM, JSW, INDLAY

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

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

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