Index: LMDZ6/trunk/libf/phylmd/radlwsw_m.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/radlwsw_m.F90	(revision 3756)
+++ LMDZ6/trunk/libf/phylmd/radlwsw_m.F90	(revision 3908)
@@ -19,6 +19,6 @@
    flag_aerosol_strat, flag_aer_feedback, &
    tau_aero, piz_aero, cg_aero,&
-   tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm,& ! rajoute par OB pour RRTM
-   tau_aero_lw_rrtm, &                                   ! rajoute par C. Kleinschmitt pour RRTM
+   tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm,& ! rajoute par OB RRTM
+   tau_aero_lw_rrtm, &              ! rajoute par C.Kleinschmitt pour RRTM
    cldtaupi, &
    qsat, flwc, fiwc, &
@@ -45,21 +45,17 @@
    ZSWFT0_i, ZFSDN0, ZFSUP0)
 
-
-
+! Modules necessaires
   USE DIMPHY
   USE assert_m, ONLY : assert
   USE infotrac_phy, ONLY : type_trac
   USE write_field_phy
+
 #ifdef REPROBUS
   USE CHEM_REP, ONLY : solaireTIME, ok_SUNTIME, ndimozon
 #endif
+
 #ifdef CPP_RRTM
 !    modules necessaires au rayonnement 
 !    -----------------------------------------
-!     USE YOMCST   , ONLY : RG       ,RD       ,RTT      ,RPI
-!     USE YOERAD   , ONLY : NSW      ,LRRTM    ,LINHOM   , LCCNL,LCCNO,
-!     USE YOERAD   , ONLY : NSW      ,LRRTM    ,LCCNL    ,LCCNO ,&
-! NSW mis dans .def MPL 20140211
-! NLW ajoute par OB
       USE YOERAD   , ONLY : NLW, LRRTM    ,LCCNL    ,LCCNO ,&
           NRADIP   , NRADLP , NICEOPT, NLIQOPT ,RCCNLND  , RCCNSEA
@@ -73,7 +69,5 @@
           RFLDD1   ,RFLDD2   ,RFLDD3   ,RFUETA   ,RASWCA,& 
           RASWCB   ,RASWCC   ,RASWCD   ,RASWCE   ,RASWCF
-!    &    RASWCB   ,RASWCC   ,RASWCD   ,RASWCE   ,RASWCF, RLINLI
       USE YOERDU   , ONLY : NUAER  ,NTRAER ,REPLOG ,REPSC  ,REPSCW ,DIFF
-!      USE YOETHF   , ONLY : RTICE
       USE YOERRTWN , ONLY : DELWAVE   ,TOTPLNK      
       USE YOMPHY3  , ONLY : RII0
@@ -81,34 +75,46 @@
       USE aero_mod
 
+! AI 02.2021
+! Besoin pour ECRAD de pctsrf, zmasq, longitude, altitude
+#ifdef CPP_ECRAD
+      USE geometry_mod, ONLY: latitude, longitude
+      USE phys_state_var_mod, ONLY: pctsrf
+      USE indice_sol_mod
+      USE time_phylmdz_mod, only: current_time
+      USE phys_cal_mod, only: day_cur
+#endif
+
   !======================================================================
   ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19960719
   ! Objet: interface entre le modele et les rayonnements
   ! Arguments:
-  ! dist-----input-R- distance astronomique terre-soleil
-  ! rmu0-----input-R- cosinus de l'angle zenithal
-  ! fract----input-R- duree d'ensoleillement normalisee
-  ! co2_ppm--input-R- concentration du gaz carbonique (en ppm)
-  ! paprs----input-R- pression a inter-couche (Pa)
-  ! pplay----input-R- pression au milieu de couche (Pa)
-  ! tsol-----input-R- temperature du sol (en K)
-  ! alb1-----input-R- albedo du sol(entre 0 et 1) dans l'interval visible 
-  ! alb2-----input-R- albedo du sol(entre 0 et 1) dans l'interval proche infra-rouge   
-  ! t--------input-R- temperature (K)
-  ! q--------input-R- vapeur d'eau (en kg/kg)
-  ! cldfra---input-R- fraction nuageuse (entre 0 et 1)
-  ! cldtaupd---input-R- epaisseur optique des nuages dans le visible (present-day value)
-  ! cldemi---input-R- emissivite des nuages dans l'IR (entre 0 et 1)
-  ! ok_ade---input-L- apply the Aerosol Direct Effect or not?
-  ! ok_aie---input-L- apply the Aerosol Indirect Effect or not?
-  ! ok_volcan-input-L- activate volcanic diags (SW heat & LW cool rate, SW & LW flux)
-  ! flag_aerosol-input-I- aerosol flag from 0 to 6
-  ! flag_aerosol_strat-input-I- use stratospheric aerosols flag (0, 1, 2)
-  ! flag_aer_feedback-input-I- activate aerosol radiative feedback (T, F)
-  ! tau_ae, piz_ae, cg_ae-input-R- aerosol optical properties (calculated in aeropt.F)
-  ! cldtaupi-input-R- epaisseur optique des nuages dans le visible
+  !                  INPUTS
+  ! dist----- input-R- distance astronomique terre-soleil
+  ! rmu0----- input-R- cosinus de l'angle zenithal
+  ! fract---- input-R- duree d'ensoleillement normalisee
+  ! co2_ppm-- input-R- concentration du gaz carbonique (en ppm)
+  ! paprs---- input-R- pression a inter-couche (Pa)
+  ! pplay---- input-R- pression au milieu de couche (Pa)
+  ! tsol----- input-R- temperature du sol (en K)
+  ! alb1----- input-R- albedo du sol(entre 0 et 1) dans l'interval visible 
+  ! alb2----- input-R- albedo du sol(entre 0 et 1) dans l'interval proche infra-rouge   
+  ! t-------- input-R- temperature (K)
+  ! q-------- input-R- vapeur d'eau (en kg/kg)
+  ! cldfra--- input-R- fraction nuageuse (entre 0 et 1)
+  ! cldtaupd- input-R- epaisseur optique des nuages dans le visible (present-day value)
+  ! cldemi--- input-R- emissivite des nuages dans l'IR (entre 0 et 1)
+  ! ok_ade--- input-L- apply the Aerosol Direct Effect or not?
+  ! ok_aie--- input-L- apply the Aerosol Indirect Effect or not?
+  ! ok_volcan input-L- activate volcanic diags (SW heat & LW cool rate, SW & LW flux)
+  ! flag_aerosol input-I- aerosol flag from 0 to 6
+  ! flag_aerosol_strat input-I- use stratospheric aerosols flag (0, 1, 2)
+  ! flag_aer_feedback  input-I- activate aerosol radiative feedback (T, F)
+  ! tau_ae, piz_ae, cg_ae input-R- aerosol optical properties (calculated in aeropt.F)
+  ! cldtaupi  input-R- epaisseur optique des nuages dans le visible
   !                   calculated for pre-industrial (pi) aerosol concentrations, i.e. with smaller
   !                   droplet concentration, thus larger droplets, thus generally cdltaupi cldtaupd
   !                   it is needed for the diagnostics of the aerosol indirect radiative forcing      
   !
+  !                  OUTPUTS
   ! heat-----output-R- echauffement atmospherique (visible) (K/jour)
   ! cool-----output-R- refroidissement dans l'IR (K/jour)
@@ -177,4 +183,8 @@
   ! 
   ! ====================================================================
+
+! ==============
+! DECLARATIONS
+! ==============
   include "YOETHF.h"
   include "YOMCST.h"
@@ -286,4 +296,41 @@
   REAL(KIND=8) PTAVE(kdlon,kflev)
   REAL(KIND=8) PWV(kdlon,kflev), PQS(kdlon,kflev)
+
+!!!!!!! Declarations specifiques pour ECRAD !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! AI 02.2021
+#ifdef CPP_ECRAD
+! ATTENTION les dimensions klon, kdlon ???
+! INPUTS
+  REAL(KIND=8) ZEMISW(klon), &              ! LW emissivity inside the window region
+               ZEMIS(klon)                  ! LW emissivity outside the window region
+  REAL(KIND=8) ZGELAM(klon), &              ! longitudes en rad
+               ZGEMU(klon)                  ! sin(latitude)
+  REAL(KIND=8) ZCO2(klon,klev), &           ! CO2 mass mixing ratios on full levels
+               ZCH4(klon,klev), &           ! CH4 mass mixing ratios on full levels
+               ZN2O(klon,klev), &           ! N2O mass mixing ratios on full levels
+               ZNO2(klon,klev), &           ! NO2 mass mixing ratios on full levels
+               ZCFC11(klon,klev), &         ! CFC11
+               ZCFC12(klon,klev), &         ! CFC12
+               ZHCFC22(klon,klev), &        ! HCFC22
+               ZCCL4(klon,klev), &          ! CCL4
+               ZO3_DP(klon,klev), ZO3_DP_i(klon,klev)            ! Ozone
+  REAL(KIND=8) ZQ_RAIN(klon,klev), &        ! Rain cloud mass mixing ratio (kg/kg) ?
+               ZQ_SNOW(klon,klev)           ! Snow cloud mass mixing ratio (kg/kg) ?
+  REAL(KIND=8) ZAEROSOL_OLD(KLON,6,KLEV), &  ! 
+               ZAEROSOL(KLON,KLEV,naero_tot) !
+! OUTPUTS
+  REAL(KIND=8) ZFLUX_DIR(klon), &           ! Direct compt of surf flux into horizontal plane
+               ZFLUX_DIR_CLEAR(klon), &     ! CS Direct 
+               ZFLUX_DIR_INTO_SUN(klon), &  ! 
+               ZFLUX_UV(klon), &            ! UV flux
+               ZFLUX_PAR(klon), &           ! photosynthetically active radiation similarly
+               ZFLUX_PAR_CLEAR(klon), &     ! CS photosynthetically 
+               ZFLUX_SW_DN_TOA(klon), &     ! DN SW flux at TOA
+               ZEMIS_OUT(klon)              ! effective broadband emissivity
+  REAL(KIND=8) ZLWDERIVATIVE(klon,klev+1)   ! LW derivatives
+  REAL(KIND=8) ZSWDIFFUSEBAND(klon,NSW), &  ! SW DN flux in diffuse albedo band 
+               ZSWDIRECTBAND(klon,NSW)      ! SW DN flux in direct albedo band
+#endif
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
   REAL(kind=8) POZON(kdlon, kflev, size(wo, 3)) ! mass fraction of ozone
@@ -401,10 +448,18 @@
   REAL zdir, zdif
 
+! =========  INITIALISATIONS ==============================================
+ IF (lldebug) THEN
+  print*,'Entree dans radlwsw '
+  print*,'************* INITIALISATIONS *****************************'
+  print*,'klon, kdlon, klev, kflev =',klon, kdlon, klev, kflev
+ ENDIF
+
   CALL assert(size(wo, 1) == klon, size(wo, 2) == klev, "radlwsw wo")
-  ! initialisation
+  
   ist=1
   iend=klon
   ktdia=1
   kmode=ist 
+! Aeros
   tauaero(:,:,:,:)=0.
   pizaero(:,:,:,:)=0.
@@ -462,23 +517,23 @@
   ENDIF
 
+ IF (lldebug) THEN
+  print*,'************** Debut boucle de 1 a ', nb_gr
+ ENDIF
+
   DO j = 1, nb_gr
     iof = kdlon*(j-1)
     DO i = 1, kdlon
       zfract(i) = fract(iof+i)
-!     zfract(i) = 1.     !!!!!!  essai MPL 19052010
       zrmu0(i) = rmu0(iof+i)
 
 
-!albedo SB >>>
-!
       IF (iflag_rrtm==0) THEN
-!
+!     Albedo
         PALBD(i,1)=alb_dif(iof+i,1)
         PALBD(i,2)=alb_dif(iof+i,2)
         PALBP(i,1)=alb_dir(iof+i,1)
         PALBP(i,2)=alb_dir(iof+i,2)
-!
-      ELSEIF (iflag_rrtm==1) THEn
-!
+! AI 02.2021 cas iflag_rrtm=1 et 2
+       ELSEIF (iflag_rrtm==1.OR.iflag_rrtm==2) THEN
         DO kk=1,NSW
           PALBD_NEW(i,kk)=alb_dif(iof+i,kk)
@@ -488,5 +543,4 @@
       ENDIF
 !albedo SB <<<
-
 
       PEMIS(i) = 1.0    !!!!! A REVOIR (MPL) 
@@ -513,4 +567,5 @@
         POZON(i,k, :) = wo(iof+i, k, :) * RG * dobson_u * 1e3 &
              / (paprs(iof+i, k) - paprs(iof+i, k+1))
+        ZO3_DP(i,k) = wo(iof+i, k, 1) * RG * dobson_u * 1e3
 !       A activer pour CCMVAL on prend l'ozone impose (MPL 07042010)
 !       POZON(i,k,:) = wo(i,k,:)  
@@ -569,4 +624,24 @@
       ENDDO
     ENDDO
+!
+! AI 02.2021
+#ifdef CPP_ECRAD
+  ZEMIS = 1.0
+  ZEMISW = 1.0
+  ZGELAM = longitude
+  ZGEMU = sin(latitude)
+  ZCO2 = RCO2 
+  ZCH4 = RCH4
+  ZN2O = RN2O
+  ZNO2 = 0.0
+  ZCFC11 = RCFC11
+  ZCFC12 = RCFC12 
+  ZHCFC22 = 0.0 
+  ZCCL4 = 0.0
+  ZQ_RAIN = 0.0
+  ZQ_SNOW = 0.0
+  ZAEROSOL_OLD = 0.0
+  ZAEROSOL = 0.0
+#endif
 !
 !===== iflag_rrtm ================================================
@@ -693,5 +768,5 @@
        ENDDO  
 !
-    ELSE  
+    ELSE IF (iflag_rrtm == 1) then  
 #ifdef CPP_RRTM
 !      if (prt_level.gt.10)write(lunout,*)'CPP_RRTM=.T.' 
@@ -804,4 +879,5 @@
             ENDDO
          ENDDO
+
 !       print *,'RADLWSW: avant RECMWFL, RI0,rmu0=',solaire,rmu0
 
@@ -850,19 +926,35 @@
 ! Nouvel appel a RECMWF (celui du cy32t0)
          CALL RECMWF_AERO (ist , iend, klon , ktdia  , klev   , kmode ,&
+!                          KST,  KEND, KPROMA, KTDIA , KLEV,  KMODE
          PALBD_NEW,PALBP_NEW, paprs_i , pplay_i , RCO2   , cldfra_i,&
+!        PALBD ,   PALBP ,    PAPRS ,   PAPRSF ,  PCCO2 , PCLFR
          POZON_i  , PAER_i  , PDP_i   , PEMIS   , rmu0   ,&
+!        PQO3  ,    PAER  ,   PDP   ,   PEMIS  , PMU0
           q_i     , qsat_i  , fiwc_i  , flwc_i  , zmasq  , t_i  ,tsol,&
+!         PQ       , PQS   , PQIWP ,     PQLWP , PSLM   , PT    , PTS, 
          ref_liq_i, ref_ice_i, &
+!        PREF_LIQ, PREF_ICE
          ref_liq_pi_i, ref_ice_pi_i, &   ! rajoute par OB pour diagnostiquer effet indirect
+!        PREF_LIQ_PI, PREF_ICE_PI
          ZEMTD_i  , ZEMTU_i , ZTRSO_i ,&
+!        PEMTD , PEMTU , PTRSO
          ZTH_i    , ZCTRSO  , ZCEMTR  , ZTRSOD  ,&
+!        PTH   ,     PCTRSO,   PCEMTR, PTRSOD
          ZLWFC    , ZLWFT_i , ZSWFC   , ZSWFT_i ,&
+!        PLWFC,     PLWFT,    PSWFC,    PSWFT, 
          PSFSWDIR , PSFSWDIF, PFSDNN  , PFSDNV  ,&
+!        PSFSWDIR,  PSFSWDIF, PFSDNN,    PFSDNV
          PPIZA_TOT, PCGA_TOT,PTAU_TOT,&
+!        PPIZA_TOT, PCGA_TOT,PTAU_TOT
          PPIZA_NAT, PCGA_NAT,PTAU_NAT,           &  ! rajoute par OB pour diagnostiquer effet direct
+!        PPIZA_NAT,PCGA_NAT,PTAU_NAT
          PTAU_LW_TOT, PTAU_LW_NAT,               &  ! rajoute par C. Kleinschmitt
+!        PTAU_LW_TOT, PTAU_LW_NAT,
          ZFLUX_i  , ZFLUC_i ,&
+!        PFLUX,      PFLUC
          ZFSDWN_i , ZFSUP_i , ZFCDWN_i, ZFCUP_i, ZFCCDWN_i, ZFCCUP_i, ZFLCCDWN_i, ZFLCCUP_i, &
+!        PFSDN ,    PFSUP ,   PFSCDN , PFSCUP,   PFSCCDN,   PFSCCUP,  PFLCCDN,    PFLCCUP
          ZTOPSWADAERO,ZSOLSWADAERO,&  ! rajoute par OB pour diagnostics
+!        PTOPSWADAERO,PSOLSWADAERO
          ZTOPSWAD0AERO,ZSOLSWAD0AERO,&
          ZTOPSWAIAERO,ZSOLSWAIAERO, &
@@ -874,4 +966,18 @@
          ZLWADAERO, & !--NL
          ok_ade, ok_aie, ok_volcan, flag_aerosol,flag_aerosol_strat, flag_aer_feedback) ! flags aerosols
+
+!--OB diagnostics
+! & PTOPSWAIAERO,PSOLSWAIAERO,&
+! & PTOPSWCFAERO,PSOLSWCFAERO,&
+! & PSWADAERO,& !--NL
+!!--LW diagnostics CK
+! & PTOPLWADAERO,PSOLLWADAERO,&
+! & PTOPLWAD0AERO,PSOLLWAD0AERO,&
+! & PTOPLWAIAERO,PSOLLWAIAERO,&
+! & PLWADAERO,& !--NL
+!!..end
+! & ok_ade, ok_aie, ok_volcan, flag_aerosol,flag_aerosol_strat,&
+! & flag_aer_feedback)
+
             
 !        print *,'RADLWSW: apres RECMWF'
@@ -1061,5 +1167,318 @@
     call abort_physic(modname, abort_message, 1)
 #endif
-    ENDIF ! iflag_rrtm
+!======================================================================
+! AI fev 2021
+    ELSE IF(iflag_rrtm == 2) THEN
+    print*,'Traitement cas iflag_rrtm = ',iflag_rrtm
+!    print*,'Mise a zero des flux '
+#ifdef CPP_ECRAD
+      DO k = 1, kflev+1
+      DO i = 1, kdlon
+        ZEMTD_i(i,k)=0.
+        ZEMTU_i(i,k)=0.
+        ZTRSO_i(i,k)=0.
+        ZTH_i(i,k)=0.
+        ZLWFT_i(i,k)=0.
+        ZSWFT_i(i,k)=0.
+        ZFLUX_i(i,1,k)=0.
+        ZFLUX_i(i,2,k)=0.
+        ZFLUC_i(i,1,k)=0.
+        ZFLUC_i(i,2,k)=0.
+        ZFSDWN_i(i,k)=0.
+        ZFCDWN_i(i,k)=0.
+        ZFCCDWN_i(i,k)=0.
+        ZFSUP_i(i,k)=0.
+        ZFCUP_i(i,k)=0.
+        ZFCCUP_i(i,k)=0.
+        ZFLCCDWN_i(i,k)=0.
+        ZFLCCUP_i(i,k)=0.
+      ENDDO
+      ENDDO
+!
+! Aerosols A REVOIR
+!      DO i = 1, kdlon
+!      DO k = 1, kflev
+!      DO kk=1, NSW
+!
+!      PTAU_TOT(i,kflev+1-k,kk)=tau_aero_sw_rrtm(i,k,2,kk)
+!      PPIZA_TOT(i,kflev+1-k,kk)=piz_aero_sw_rrtm(i,k,2,kk)
+!      PCGA_TOT(i,kflev+1-k,kk)=cg_aero_sw_rrtm(i,k,2,kk)
+!
+!      PTAU_NAT(i,kflev+1-k,kk)=tau_aero_sw_rrtm(i,k,1,kk)
+!      PPIZA_NAT(i,kflev+1-k,kk)=piz_aero_sw_rrtm(i,k,1,kk)
+!      PCGA_NAT(i,kflev+1-k,kk)=cg_aero_sw_rrtm(i,k,1,kk)
+!
+!      ENDDO
+!      ENDDO
+!      ENDDO
+!-end OB
+!
+!      DO i = 1, kdlon
+!      DO k = 1, kflev
+!      DO kk=1, NLW
+!
+!      PTAU_LW_TOT(i,kflev+1-k,kk)=tau_aero_lw_rrtm(i,k,2,kk)
+!      PTAU_LW_NAT(i,kflev+1-k,kk)=tau_aero_lw_rrtm(i,k,1,kk)
+!
+!      ENDDO
+!      ENDDO
+!      ENDDO
+!-end C. Kleinschmitt
+!      
+      DO i = 1, kdlon
+      ZCTRSO(i,1)=0.
+      ZCTRSO(i,2)=0.
+      ZCEMTR(i,1)=0.
+      ZCEMTR(i,2)=0.
+      ZTRSOD(i)=0.
+      ZLWFC(i,1)=0.
+      ZLWFC(i,2)=0.
+      ZSWFC(i,1)=0.
+      ZSWFC(i,2)=0.
+      PFSDNN(i)=0.
+      PFSDNV(i)=0.
+      DO kk = 1, NSW
+        PSFSWDIR(i,kk)=0.
+        PSFSWDIF(i,kk)=0.
+      ENDDO
+      ENDDO
+!----- Fin des mises a zero des tableaux output -------------------              
+
+! On met les donnees dans l'ordre des niveaux ecrad
+!         print*,'On inverse sur la verticale '
+         paprs_i(:,1)=paprs(:,klev+1)
+         DO k=1,klev
+            paprs_i(1:klon,k+1) =paprs(1:klon,klev+1-k)
+            pplay_i(1:klon,k)   =pplay(1:klon,klev+1-k)
+            cldfra_i(1:klon,k)  =cldfra(1:klon,klev+1-k)
+            PDP_i(1:klon,k)     =PDP(1:klon,klev+1-k)
+            t_i(1:klon,k)       =t(1:klon,klev+1-k)
+            q_i(1:klon,k)       =q(1:klon,klev+1-k)
+            qsat_i(1:klon,k)    =qsat(1:klon,klev+1-k)
+            flwc_i(1:klon,k)    =flwc(1:klon,klev+1-k)
+            fiwc_i(1:klon,k)    =fiwc(1:klon,klev+1-k)
+            ref_liq_i(1:klon,k) =ref_liq(1:klon,klev+1-k)
+            ref_ice_i(1:klon,k) =ref_ice(1:klon,klev+1-k)
+!-OB
+            ref_liq_pi_i(1:klon,k) =ref_liq_pi(1:klon,klev+1-k)
+            ref_ice_pi_i(1:klon,k) =ref_ice_pi(1:klon,klev+1-k)
+         ENDDO
+         DO k=1,kflev
+            ZO3_DP_i(1:klon,k)=ZO3_DP(1:klon,kflev+1-k)
+            DO i=1,6
+            PAER_i(1:klon,k,i)=PAER(1:klon,kflev+1-k,i)
+            ENDDO
+         ENDDO
+! AI 02.2021
+! Calcul of ZTH_i
+         DO K=2,KLEV
+            ZTH_i(:,K)=&
+              & (t_i(:,K-1)*pplay_i(:,K-1)*(pplay_i(:,K)-paprs_i(:,K))&
+              & +t_i(:,K)*pplay_i(:,K)*(paprs_i(:,K)-pplay_i(:,K-1)))&
+              & *(1.0/(paprs_i(:,K)*(pplay_i(:,K)-pplay_i(:,K-1))))
+         ENDDO
+            ZTH_i(:,KLEV+1)=tsol(:)
+            ZTH_i(:,1)=t_i(:,1)-pplay_i(:,1)*(t_i(:,1)-ZTH_i(:,2))&
+                      & /(pplay_i(:,1)-paprs_i(:,2))
+
+      print *,'RADLWSW: avant RADIATION_SCHEME '
+      IF (lldebug) THEN
+        CALL writefield_phy('rmu0',rmu0,1)
+        CALL writefield_phy('tsol',tsol,1)
+        CALL writefield_phy('emissiv_out',ZEMIS,1)
+        CALL writefield_phy('emissiv_in',ZEMISW,1)
+        CALL writefield_phy('pctsrf_ter',pctsrf(:,is_ter),1)
+        CALL writefield_phy('pctsrf_oce',pctsrf(:,is_oce),1)
+        CALL writefield_phy('ZGELAM',ZGELAM,1)
+        CALL writefield_phy('ZGEMU',ZGEMU,1)
+        CALL writefield_phy('zmasq',zmasq,1)
+        CALL writefield_phy('paprs_i',paprs_i,klev+1)
+        CALL writefield_phy('pplay_i',pplay_i,klev)
+        CALL writefield_phy('t_i',t_i,klev)
+        CALL writefield_phy('ZTH_i',ZTH_i,klev+1)
+        CALL writefield_phy('cldfra_i',cldfra_i,klev)
+        CALL writefield_phy('paer_i',PAER_i,klev)
+        CALL writefield_phy('q_i',q_i,klev)
+        CALL writefield_phy('fiwc_i',fiwc_i,klev)
+        CALL writefield_phy('flwc_i',flwc_i,klev)
+        CALL writefield_phy('palbd_new',PALBD_NEW,NSW)
+        CALL writefield_phy('palbp_new',PALBP_NEW,NSW)
+        CALL writefield_phy('ZO3_DP',ZO3_DP,klev)
+      ENDIF
+
+      CALL RADIATION_SCHEME &
+      & (ist, iend, klon, klev, naero_tot, NSW, &
+! ??? naero_tot
+      & day_cur, current_time, & 
+      & solaire, &
+      & rmu0, tsol, PALBD_NEW,PALBP_NEW, &   
+!       PEMIS_WINDOW (???), &
+      &  ZEMIS, ZEMISW, &
+!       PCCN_LAND, PCCN_SEA, & ???
+      & pctsrf(:,is_ter), pctsrf(:,is_oce), &
+!       longitude(rad), sin(latitude), PMASQ_ ???
+      & ZGELAM, ZGEMU, zmasq, &
+      & pplay_i, t_i, &
+!       PTEMPERATURE_H ?, &
+      & paprs_i, ZTH_i, q_i, qsat_i, & 
+      & ZCO2, ZCH4, ZN2O, ZNO2, ZCFC11, ZCFC12, ZHCFC22, ZCCL4, ZO3_DP_i, &
+      & cldfra_i, flwc_i, fiwc_i, ZQ_RAIN, ZQ_SNOW, &  
+      & ref_liq_i, ref_ice_i, &
+!       aerosols
+      & ZAEROSOL_OLD, ZAEROSOL, &
+! Outputs
+!       Net flux
+      & ZSWFT_i, ZLWFT_i, ZSWFT0_i, ZLWFT0_i, &
+!       DN flux
+      & ZFSDWN_i, ZFLUX_i(:,2,:), ZFCDWN_i, ZFLUC_i(:,2,:), &
+!       UP flux
+      & ZFSUP_i, ZFLUX_i(:,1,:), ZFCUP_i, ZFLUC_i(:,1,:), &
+!       Surf Direct flux
+      & ZFLUX_DIR, ZFLUX_DIR_CLEAR, ZFLUX_DIR_INTO_SUN, &
+!       UV and para flux
+      & ZFLUX_UV, ZFLUX_PAR, ZFLUX_PAR_CLEAR, &
+!      & ZFLUX_SW_DN_TOA, 
+      & ZEMIS_OUT, ZLWDERIVATIVE, &
+      & PSFSWDIF, PSFSWDIR)
+
+      print *,'========= RADLWSW: apres RADIATION_SCHEME ==================== '
+  IF (lldebug) THEN
+    if (0.eq.0) then
+      print *,' Net Flux '
+      print *,'ZSWFT_i =', ZSWFT_i
+      print *,'ZLWFT_i =', ZLWFT_i
+      print *,'ZSWFT0_i =', ZSWFT0_i
+      print*,'ZLWFT0_i =', ZLWFT0_i
+
+      print*,'DN Flux '
+      print*,'ZFSDWN_i =', ZFSDWN_i
+      print*,'ZFLUX_i(:,2,:)', ZFLUX_i(:,2,:)
+      print*,'ZFCDWN_i =', ZFCDWN_i
+      print*,'ZFLUC_i(:,2,:) =', ZFLUC_i(:,2,:) 
+
+      print*,'UP Flux '
+      print*,'ZFSUP_i =', ZFSUP_i 
+      print*,'ZFLUX_i(:,1,:) =', ZFLUX_i(:,1,:)
+      print*,'ZFCUP_i =', ZFCUP_i
+      print*,'ZFLUC_i(:,1,:) =', ZFLUC_i(:,1,:)
+
+      print*,'UV and para flux '
+      print*,'ZFLUX_DIR =', ZFLUX_DIR
+      print*,'ZFLUX_DIR_CLEAR', ZFLUX_DIR_CLEAR
+      print*,'ZFLUX_DIR_INTO_SUN', ZFLUX_DIR_INTO_SUN
+    endif
+   ENDIF
+
+      IF (lldebug) THEN
+        CALL writefield_phy('zlwft_i',ZLWFT_i,klev+1)
+        CALL writefield_phy('zlwft0_i',ZLWFT0_i,klev+1)
+        CALL writefield_phy('zswft_i',ZSWFT_i,klev+1)
+        CALL writefield_phy('zswft0_i',ZSWFT0_i,klev+1)
+        CALL writefield_phy('psfswdir',PSFSWDIR,6)
+        CALL writefield_phy('psfswdif',PSFSWDIF,6)
+        CALL writefield_phy('zflux_i',ZFLUX_i,klev+1)
+        CALL writefield_phy('zfluc_i',ZFLUC_i,klev+1)
+        CALL writefield_phy('zfsdwn_i',ZFSDWN_i,klev+1)
+        CALL writefield_phy('zfsup_i',ZFSUP_i,klev+1)
+        CALL writefield_phy('zfcdwn_i',ZFCDWN_i,klev+1)
+        CALL writefield_phy('zfcup_i',ZFCUP_i,klev+1)
+      ENDIF
+! ---------
+! On retablit l'ordre des niveaux lmd pour les tableaux de sortie
+! D autre part, on multiplie les resultats SW par fract pour etre coherent
+! avec l ancien rayonnement AR4. Si nuit, fract=0 donc pas de 
+! rayonnement SW. (MPL 260609)
+      print*,'On retablit ordre des niveaux verticaux'
+      print*,'On multiplie les flux SW par fract ?'
+      DO k=0,klev
+         DO i=1,klon
+         ZEMTD(i,k+1)  = ZEMTD_i(i,k+1)
+         ZEMTU(i,k+1)  = ZEMTU_i(i,k+1)
+         ZTRSO(i,k+1)  = ZTRSO_i(i,k+1)
+         ZTH(i,k+1)    = ZTH_i(i,k+1)
+         ZFLUP(i,k+1)  = ZFLUX_i(i,1,k+1)
+         ZFLDN(i,k+1)  = ZFLUX_i(i,2,k+1)
+         ZFLUP0(i,k+1) = ZFLUC_i(i,1,k+1)
+         ZFLDN0(i,k+1) = ZFLUC_i(i,2,k+1)
+         ZFSDN(i,k+1)  = ZFSDWN_i(i,k+1)*fract(i)
+         ZFSDN0(i,k+1) = ZFCDWN_i(i,k+1)*fract(i)
+         ZFSDNC0(i,k+1)= ZFCCDWN_i(i,k+1)*fract(i)
+         ZFSUP (i,k+1) = ZFSUP_i(i,k+1)*fract(i)
+         ZFSUP0(i,k+1) = ZFCUP_i(i,k+1)*fract(i)
+         ZFSUPC0(i,k+1)= ZFCCUP_i(i,k+1)*fract(i)
+         ZFLDNC0(i,k+1)= ZFLCCDWN_i(i,k+1)
+         ZFLUPC0(i,k+1)= ZFLCCUP_i(i,k+1)
+         IF (ok_volcan) THEN
+            ZSWADAERO(i,k+1)=ZSWADAERO(i,k+1)*fract(i) !--NL
+         ENDIF
+         
+!   Nouveau calcul car visiblement ZSWFT et ZSWFC sont nuls dans RRTM cy32
+!   en sortie de radlsw.F90 - MPL 7.01.09
+         ZSWFT(i,k+1)  = (ZFSDWN_i(i,k+1)-ZFSUP_i(i,k+1))*fract(i)
+         ZSWFT0_i(i,k+1) = (ZFCDWN_i(i,k+1)-ZFCUP_i(i,k+1))*fract(i)
+         ZLWFT(i,k+1) =-ZFLUX_i(i,2,k+1)-ZFLUX_i(i,1,k+1)
+         ZLWFT0_i(i,k+1)=-ZFLUC_i(i,2,k+1)-ZFLUC_i(i,1,k+1)
+         ENDDO
+      ENDDO
+
+!--ajout OB
+      ZTOPSWADAERO(:) =ZTOPSWADAERO(:) *fract(:)
+      ZSOLSWADAERO(:) =ZSOLSWADAERO(:) *fract(:)
+      ZTOPSWAD0AERO(:)=ZTOPSWAD0AERO(:)*fract(:)
+      ZSOLSWAD0AERO(:)=ZSOLSWAD0AERO(:)*fract(:)
+      ZTOPSWAIAERO(:) =ZTOPSWAIAERO(:) *fract(:)
+      ZSOLSWAIAERO(:) =ZSOLSWAIAERO(:) *fract(:)
+      ZTOPSWCF_AERO(:,1)=ZTOPSWCF_AERO(:,1)*fract(:) 
+      ZTOPSWCF_AERO(:,2)=ZTOPSWCF_AERO(:,2)*fract(:) 
+      ZTOPSWCF_AERO(:,3)=ZTOPSWCF_AERO(:,3)*fract(:) 
+      ZSOLSWCF_AERO(:,1)=ZSOLSWCF_AERO(:,1)*fract(:)
+      ZSOLSWCF_AERO(:,2)=ZSOLSWCF_AERO(:,2)*fract(:)
+      ZSOLSWCF_AERO(:,3)=ZSOLSWCF_AERO(:,3)*fract(:)
+
+! ---------
+! On renseigne les champs LMDz, pour avoir la meme chose qu'en sortie de
+! LW_LMDAR4 et SW_LMDAR4
+
+      !--fraction of diffuse radiation in surface SW downward radiation
+      DO i = 1, kdlon
+       IF (fract(i).GT.0.0) THEN
+         zdir=SUM(PSFSWDIR(i,:))
+         zdif=SUM(PSFSWDIF(i,:))
+         zsolswfdiff(i) = zdif/(zdir+zdif)
+       ELSE  !--night
+         zsolswfdiff(i) = 1.0
+       ENDIF
+      ENDDO 
+!
+      DO i = 1, kdlon
+         zsolsw(i)    = ZSWFT(i,1)
+         zsolsw0(i)   = ZSWFT0_i(i,1)
+         ztopsw(i)    = ZSWFT(i,klev+1)
+         ztopsw0(i)   = ZSWFT0_i(i,klev+1)
+         zsollw(i)    = ZLWFT(i,1)
+         zsollw0(i)   = ZLWFT0_i(i,1)
+         ztoplw(i)    = ZLWFT(i,klev+1)*(-1)
+         ztoplw0(i)   = ZLWFT0_i(i,klev+1)*(-1)
+!         
+         zsollwdown(i)= -1.*ZFLDN(i,1)
+      ENDDO
+
+      DO k=1,kflev
+         DO i=1,kdlon
+           zheat(i,k)=(ZSWFT(i,k+1)-ZSWFT(i,k))*RDAY*RG/RCPD/PDP(i,k)
+           zheat0(i,k)=(ZSWFT0_i(i,k+1)-ZSWFT0_i(i,k))*RDAY*RG/RCPD/PDP(i,k)
+           zcool(i,k)=(ZLWFT(i,k)-ZLWFT(i,k+1))*RDAY*RG/RCPD/PDP(i,k)
+           zcool0(i,k)=(ZLWFT0_i(i,k)-ZLWFT0_i(i,k+1))*RDAY*RG/RCPD/PDP(i,k)
+           IF (ok_volcan) THEN
+              zheat_volc(i,k)=(ZSWADAERO(i,k+1)-ZSWADAERO(i,k))*RG/RCPD/PDP(i,k) !NL
+              zcool_volc(i,k)=(ZLWADAERO(i,k)-ZLWADAERO(i,k+1))*RG/RCPD/PDP(i,k) !NL
+           ENDIF
+         ENDDO
+      ENDDO
+#endif  
+  print*,'Fin traitement ECRAD'
+! Fin ECRAD
+  ENDIF        ! iflag_rrtm
+! ecrad
 !======================================================================
 
@@ -1102,9 +1521,4 @@
           solswad_aero(iof+i) = zsolswadaero(i)
           solswad0_aero(iof+i) = zsolswad0aero(i)
-! MS the following lines seem to be wrong, why is iof on right hand side???
-!          topsw_aero(iof+i,:) = ztopsw_aero(iof+i,:)
-!          topsw0_aero(iof+i,:) = ztopsw0_aero(iof+i,:)
-!          solsw_aero(iof+i,:) = zsolsw_aero(iof+i,:)
-!          solsw0_aero(iof+i,:) = zsolsw0_aero(iof+i,:)
           topsw_aero(iof+i,:) = ztopsw_aero(i,:)
           topsw0_aero(iof+i,:) = ztopsw0_aero(i,:)
