Changeset 4866


Ignore:
Timestamp:
Mar 21, 2024, 10:03:50 PM (7 weeks ago)
Author:
lguez
Message:

Indent file

File:
1 edited

Legend:

Unmodified
Added
Removed
  • LMDZ6/trunk/libf/phylmd/radlwsw_m.F90

    r4864 r4866  
    88contains
    99
    10 SUBROUTINE radlwsw( &
    11    debut, dist, rmu0, fract, &
    12 !albedo SB >>>
    13 !  paprs, pplay,tsol,alb1, alb2, &
    14    paprs, pplay,tsol,SFRWL,alb_dir, alb_dif, &
    15 !albedo SB <<<
    16    t,q,wo,&
    17    cldfra, cldemi, cldtaupd,&
    18    ok_ade, ok_aie, ok_volcan, flag_volc_surfstrat, flag_aerosol,&
    19    flag_aerosol_strat, flag_aer_feedback, &
    20    tau_aero, piz_aero, cg_aero,&
    21    tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm,& ! rajoute par OB RRTM
    22    tau_aero_lw_rrtm, &              ! rajoute par C.Kleinschmitt pour RRTM
    23    cldtaupi, &
    24    qsat, flwc, fiwc, &
    25    ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
    26    namelist_ecrad_file, &
    27    heat,heat0,cool,cool0,albpla,&
    28    heat_volc, cool_volc,&
    29    topsw,toplw,solsw,solswfdiff,sollw,&
    30    sollwdown,&
    31    topsw0,toplw0,solsw0,sollw0,&
    32    lwdnc0, lwdn0, lwdn, lwupc0, lwup0, lwup,&
    33    swdnc0, swdn0, swdn, swupc0, swup0, swup,&
    34    topswad_aero, solswad_aero,&
    35    topswai_aero, solswai_aero, &
    36    topswad0_aero, solswad0_aero,&
    37    topsw_aero, topsw0_aero,&
    38    solsw_aero, solsw0_aero, &
    39    topswcf_aero, solswcf_aero,&
    40 !-C. Kleinschmitt for LW diagnostics
    41    toplwad_aero, sollwad_aero,&
    42    toplwai_aero, sollwai_aero, &
    43    toplwad0_aero, sollwad0_aero, &
    44 !-end
    45    ZLWFT0_i, ZFLDN0, ZFLUP0, &
    46    ZSWFT0_i, ZFSDN0, ZFSUP0, &
    47    cloud_cover_sw)
    48 
    49 ! Modules necessaires
    50   USE DIMPHY
    51   USE assert_m, ONLY : assert
    52   USE infotrac_phy, ONLY : type_trac
    53   USE write_field_phy
     10  SUBROUTINE radlwsw( &
     11       debut, dist, rmu0, fract, &
     12       !albedo SB >>>
     13       !  paprs, pplay,tsol,alb1, alb2, &
     14       paprs, pplay,tsol,SFRWL,alb_dir, alb_dif, &
     15       !albedo SB <<<
     16       t,q,wo,&
     17       cldfra, cldemi, cldtaupd,&
     18       ok_ade, ok_aie, ok_volcan, flag_volc_surfstrat, flag_aerosol,&
     19       flag_aerosol_strat, flag_aer_feedback, &
     20       tau_aero, piz_aero, cg_aero,&
     21       tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm,& ! rajoute par OB RRTM
     22       tau_aero_lw_rrtm, &              ! rajoute par C.Kleinschmitt pour RRTM
     23       cldtaupi, &
     24       qsat, flwc, fiwc, &
     25       ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
     26       namelist_ecrad_file, &
     27       heat,heat0,cool,cool0,albpla,&
     28       heat_volc, cool_volc,&
     29       topsw,toplw,solsw,solswfdiff,sollw,&
     30       sollwdown,&
     31       topsw0,toplw0,solsw0,sollw0,&
     32       lwdnc0, lwdn0, lwdn, lwupc0, lwup0, lwup,&
     33       swdnc0, swdn0, swdn, swupc0, swup0, swup,&
     34       topswad_aero, solswad_aero,&
     35       topswai_aero, solswai_aero, &
     36       topswad0_aero, solswad0_aero,&
     37       topsw_aero, topsw0_aero,&
     38       solsw_aero, solsw0_aero, &
     39       topswcf_aero, solswcf_aero,&
     40       !-C. Kleinschmitt for LW diagnostics
     41       toplwad_aero, sollwad_aero,&
     42       toplwai_aero, sollwai_aero, &
     43       toplwad0_aero, sollwad0_aero, &
     44       !-end
     45       ZLWFT0_i, ZFLDN0, ZFLUP0, &
     46       ZSWFT0_i, ZFSDN0, ZFSUP0, &
     47       cloud_cover_sw)
     48
     49    ! Modules necessaires
     50    USE DIMPHY
     51    USE assert_m, ONLY : assert
     52    USE infotrac_phy, ONLY : type_trac
     53    USE write_field_phy
    5454
    5555#ifdef REPROBUS
    56   USE CHEM_REP, ONLY : solaireTIME, ok_SUNTIME, ndimozon
     56    USE CHEM_REP, ONLY : solaireTIME, ok_SUNTIME, ndimozon
    5757#endif
    5858
    5959#ifdef CPP_RRTM
    60 !    modules necessaires au rayonnement
    61 !    -----------------------------------------
    62       USE YOERAD   , ONLY : NLW, LRRTM    ,LCCNL    ,LCCNO ,&
    63           NRADIP   , NRADLP , NICEOPT, NLIQOPT ,RCCNLND  , RCCNSEA
    64       USE YOELW    , ONLY : NSIL     ,NTRA     ,NUA      ,TSTAND   ,XP
    65       USE YOESW    , ONLY : RYFWCA   ,RYFWCB   ,RYFWCC   ,RYFWCD,&   
    66           RYFWCE   ,RYFWCF   ,REBCUA   ,REBCUB   ,REBCUC,&   
    67           REBCUD   ,REBCUE   ,REBCUF   ,REBCUI   ,REBCUJ,& 
    68           REBCUG   ,REBCUH   ,RHSAVI   ,RFULIO   ,RFLAA0,& 
    69           RFLAA1   ,RFLBB0   ,RFLBB1   ,RFLBB2   ,RFLBB3,& 
    70           RFLCC0   ,RFLCC1   ,RFLCC2   ,RFLCC3   ,RFLDD0,& 
    71           RFLDD1   ,RFLDD2   ,RFLDD3   ,RFUETA   ,RASWCA,&
    72           RASWCB   ,RASWCC   ,RASWCD   ,RASWCE   ,RASWCF
    73       USE YOERDU   , ONLY : NUAER  ,NTRAER ,REPLOG ,REPSC  ,REPSCW ,DIFF
    74       USE YOERRTWN , ONLY : DELWAVE   ,TOTPLNK     
    75       USE YOMPHY3  , ONLY : RII0
     60    !    modules necessaires au rayonnement
     61    !    -----------------------------------------
     62    USE YOERAD   , ONLY : NLW, LRRTM    ,LCCNL    ,LCCNO ,&
     63         NRADIP   , NRADLP , NICEOPT, NLIQOPT ,RCCNLND  , RCCNSEA
     64    USE YOELW    , ONLY : NSIL     ,NTRA     ,NUA      ,TSTAND   ,XP
     65    USE YOESW    , ONLY : RYFWCA   ,RYFWCB   ,RYFWCC   ,RYFWCD,&   
     66         RYFWCE   ,RYFWCF   ,REBCUA   ,REBCUB   ,REBCUC,&   
     67         REBCUD   ,REBCUE   ,REBCUF   ,REBCUI   ,REBCUJ,& 
     68         REBCUG   ,REBCUH   ,RHSAVI   ,RFULIO   ,RFLAA0,& 
     69         RFLAA1   ,RFLBB0   ,RFLBB1   ,RFLBB2   ,RFLBB3,& 
     70         RFLCC0   ,RFLCC1   ,RFLCC2   ,RFLCC3   ,RFLDD0,& 
     71         RFLDD1   ,RFLDD2   ,RFLDD3   ,RFUETA   ,RASWCA,&
     72         RASWCB   ,RASWCC   ,RASWCD   ,RASWCE   ,RASWCF
     73    USE YOERDU   , ONLY : NUAER  ,NTRAER ,REPLOG ,REPSC  ,REPSCW ,DIFF
     74    USE YOERRTWN , ONLY : DELWAVE   ,TOTPLNK     
     75    USE YOMPHY3  , ONLY : RII0
    7676#endif
    77       USE aero_mod
    78 
    79 ! AI 02.2021
    80 ! Besoin pour ECRAD de pctsrf, zmasq, longitude, altitude
     77    USE aero_mod
     78
     79    ! AI 02.2021
     80    ! Besoin pour ECRAD de pctsrf, zmasq, longitude, altitude
    8181#ifdef CPP_ECRAD
    82       USE phys_local_var_mod, ONLY: rhcl, m_allaer
    83       USE geometry_mod, ONLY: latitude, longitude
    84       USE phys_state_var_mod, ONLY: pctsrf
    85       USE indice_sol_mod
    86       USE time_phylmdz_mod, only: current_time
    87       USE phys_cal_mod, only: day_cur
    88       USE interface_lmdz_ecrad
     82    USE phys_local_var_mod, ONLY: rhcl, m_allaer
     83    USE geometry_mod, ONLY: latitude, longitude
     84    USE phys_state_var_mod, ONLY: pctsrf
     85    USE indice_sol_mod
     86    USE time_phylmdz_mod, only: current_time
     87    USE phys_cal_mod, only: day_cur
     88    USE interface_lmdz_ecrad
    8989#endif
    9090
    91   !======================================================================
    92   ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19960719
    93   ! Objet: interface entre le modele et les rayonnements
    94   ! Arguments:
    95   !                  INPUTS
    96   ! dist----- input-R- distance astronomique terre-soleil
    97   ! rmu0----- input-R- cosinus de l'angle zenithal
    98   ! fract---- input-R- duree d'ensoleillement normalisee
    99   ! co2_ppm-- input-R- concentration du gaz carbonique (en ppm)
    100   ! paprs---- input-R- pression a inter-couche (Pa)
    101   ! pplay---- input-R- pression au milieu de couche (Pa)
    102   ! tsol----- input-R- temperature du sol (en K)
    103   ! alb1----- input-R- albedo du sol(entre 0 et 1) dans l'interval visible
    104   ! alb2----- input-R- albedo du sol(entre 0 et 1) dans l'interval proche infra-rouge   
    105   ! t-------- input-R- temperature (K)
    106   ! q-------- input-R- vapeur d'eau (en kg/kg)
    107   ! cldfra--- input-R- fraction nuageuse (entre 0 et 1)
    108   ! cldtaupd- input-R- epaisseur optique des nuages dans le visible (present-day value)
    109   ! cldemi--- input-R- emissivite des nuages dans l'IR (entre 0 et 1)
    110   ! ok_ade--- input-L- apply the Aerosol Direct Effect or not?
    111   ! ok_aie--- input-L- apply the Aerosol Indirect Effect or not?
    112   ! ok_volcan input-L- activate volcanic diags (SW heat & LW cool rate, SW & LW flux)
    113   ! flag_volc_surfstrat input-I- activate volcanic surf cooling or strato heating (or nothing)
    114   ! flag_aerosol input-I- aerosol flag from 0 to 6
    115   ! flag_aerosol_strat input-I- use stratospheric aerosols flag (0, 1, 2)
    116   ! flag_aer_feedback  input-I- activate aerosol radiative feedback (T, F)
    117   ! tau_ae, piz_ae, cg_ae input-R- aerosol optical properties (calculated in aeropt.F)
    118   ! cldtaupi  input-R- epaisseur optique des nuages dans le visible
    119   !                   calculated for pre-industrial (pi) aerosol concentrations, i.e. with smaller
    120   !                   droplet concentration, thus larger droplets, thus generally cdltaupi cldtaupd
    121   !                   it is needed for the diagnostics of the aerosol indirect radiative forcing     
    122   !
    123   !                  OUTPUTS
    124   ! heat-----output-R- echauffement atmospherique (visible) (K/jour)
    125   ! cool-----output-R- refroidissement dans l'IR (K/jour)
    126   ! albpla---output-R- albedo planetaire (entre 0 et 1)
    127   ! topsw----output-R- flux solaire net au sommet de l'atm.
    128   ! toplw----output-R- ray. IR montant au sommet de l'atmosphere
    129   ! solsw----output-R- flux solaire net a la surface
    130   ! solswfdiff----output-R- fraction de rayonnement diffus pour le flux solaire descendant a la surface
    131   ! sollw----output-R- ray. IR montant a la surface
    132   ! solswad---output-R- ray. solaire net absorbe a la surface (aerosol dir)
    133   ! topswad---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol dir)
    134   ! solswai---output-R- ray. solaire net absorbe a la surface (aerosol ind)
    135   ! topswai---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol ind)
    136   !
    137   ! heat_volc-----output-R- echauffement atmospherique  du au forcage volcanique (visible) (K/s)
    138   ! cool_volc-----output-R- refroidissement dans l'IR du au forcage volcanique (K/s)
    139   !
    140   ! ATTENTION: swai and swad have to be interpreted in the following manner:
    141   ! ---------
    142   ! ok_ade=F & ok_aie=F -both are zero
    143   ! ok_ade=T & ok_aie=F -aerosol direct forcing is F_{AD} = topsw-topswad
    144   !                        indirect is zero
    145   ! ok_ade=F & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai
    146   !                        direct is zero
    147   ! ok_ade=T & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai
    148   !                        aerosol direct forcing is F_{AD} = topswai-topswad
    149   !
    150   ! --------- RRTM: output RECMWFL
    151   ! ZEMTD (KPROMA,KLEV+1)         ; TOTAL DOWNWARD LONGWAVE EMISSIVITY
    152   ! ZEMTU (KPROMA,KLEV+1)         ; TOTAL UPWARD   LONGWAVE EMISSIVITY
    153   ! ZTRSO (KPROMA,KLEV+1)         ; TOTAL SHORTWAVE TRANSMISSIVITY
    154   ! ZTH   (KPROMA,KLEV+1)         ; HALF LEVEL TEMPERATURE
    155   ! ZCTRSO(KPROMA,2)              ; CLEAR-SKY SHORTWAVE TRANSMISSIVITY
    156   ! ZCEMTR(KPROMA,2)              ; CLEAR-SKY NET LONGWAVE EMISSIVITY
    157   ! ZTRSOD(KPROMA)                ; TOTAL-SKY SURFACE SW TRANSMISSITY
    158   ! ZLWFC (KPROMA,2)              ; CLEAR-SKY LONGWAVE FLUXES
    159   ! ZLWFT (KPROMA,KLEV+1)         ; TOTAL-SKY LONGWAVE FLUXES
    160   ! ZLWFT0(KPROMA,KLEV+1)         ; CLEAR-SKY LONGWAVE FLUXES      ! added by MPL 090109
    161   ! ZSWFC (KPROMA,2)              ; CLEAR-SKY SHORTWAVE FLUXES
    162   ! ZSWFT (KPROMA,KLEV+1)         ; TOTAL-SKY SHORTWAVE FLUXES
    163   ! ZSWFT0(KPROMA,KLEV+1)         ; CLEAR-SKY SHORTWAVE FLUXES     ! added by MPL 090109
    164   ! ZFLUX (KLON,2,KLEV+1)         ; TOTAL LW FLUXES  1=up, 2=DWN   ! added by MPL 080411
    165   ! ZFLUC (KLON,2,KLEV+1)         ; CLEAR SKY LW FLUXES            ! added by MPL 080411
    166   ! ZFSDWN(klon,KLEV+1)           ; TOTAL SW  DWN FLUXES           ! added by MPL 080411
    167   ! ZFCDWN(klon,KLEV+1)           ; CLEAR SKY SW  DWN FLUXES       ! added by MPL 080411
    168   ! ZFCCDWN(klon,KLEV+1)          ; CLEAR SKY CLEAN (NO AEROSOL) SW  DWN FLUXES      ! added by OB 211117
    169   ! ZFSUP (klon,KLEV+1)           ; TOTAL SW  UP  FLUXES           ! added by MPL 080411
    170   ! ZFCUP (klon,KLEV+1)           ; CLEAR SKY SW  UP  FLUXES       ! added by MPL 080411
    171   ! ZFCCUP (klon,KLEV+1)          ; CLEAR SKY CLEAN (NO AEROSOL) SW  UP  FLUXES      ! added by OB 211117
    172   ! ZFLCCDWN(klon,KLEV+1)         ; CLEAR SKY CLEAN (NO AEROSOL) LW  DWN FLUXES      ! added by OB 211117
    173   ! ZFLCCUP (klon,KLEV+1)         ; CLEAR SKY CLEAN (NO AEROSOL) LW  UP  FLUXES      ! added by OB 211117
    174  
    175   !======================================================================
    176  
    177   ! ====================================================================
    178   ! Adapte au modele de chimie INCA par Celine Deandreis & Anne Cozic -- 2009
    179   ! 1 = ZERO   
    180   ! 2 = AER total   
    181   ! 3 = NAT   
    182   ! 4 = BC   
    183   ! 5 = SO4   
    184   ! 6 = POM   
    185   ! 7 = DUST   
    186   ! 8 = SS   
    187   ! 9 = NO3   
    188   !
    189   ! ====================================================================
    190 
    191 ! ==============
    192 ! DECLARATIONS
    193 ! ==============
    194   include "YOETHF.h"
    195   include "YOMCST.h"
    196   include "clesphys.h"
    197 
    198 ! Input arguments
    199   REAL,    INTENT(in)  :: dist
    200   REAL,    INTENT(in)  :: rmu0(KLON), fract(KLON)
    201   REAL,    INTENT(in)  :: paprs(KLON,KLEV+1), pplay(KLON,KLEV)
    202 !albedo SB >>>
    203 ! REAL,    INTENT(in)  :: alb1(KLON), alb2(KLON), tsol(KLON)
    204   REAL,    INTENT(in)  :: tsol(KLON)
    205   REAL,    INTENT(in) :: alb_dir(KLON,NSW),alb_dif(KLON,NSW)
    206   REAL,    INTENT(in) :: SFRWL(6)
    207 !albedo SB <<<
    208   REAL,    INTENT(in)  :: t(KLON,KLEV), q(KLON,KLEV)
    209 
    210   REAL, INTENT(in):: wo(:, :, :) ! dimension(KLON,KLEV, 1 or 2)
    211   ! column-density of ozone in a layer, in kilo-Dobsons
    212   ! "wo(:, :, 1)" is for the average day-night field,
    213   ! "wo(:, :, 2)" is for daylight time.
    214 
    215   LOGICAL, INTENT(in)  :: ok_ade, ok_aie                                 ! switches whether to use aerosol direct (indirect) effects or not
    216   LOGICAL, INTENT(in)  :: ok_volcan                                      ! produce volcanic diags (SW/LW heat flux and rate)
    217   INTEGER, INTENT(in)  :: flag_volc_surfstrat                            ! allow to impose volcanic cooling rate at surf or heating in strato
    218   LOGICAL              :: lldebug=.false.
    219   INTEGER, INTENT(in)  :: flag_aerosol                                   ! takes value 0 (no aerosol) or 1 to 6 (aerosols)
    220   INTEGER, INTENT(in)  :: flag_aerosol_strat                             ! use stratospheric aerosols
    221   LOGICAL, INTENT(in)  :: flag_aer_feedback                              ! activate aerosol radiative feedback
    222   REAL,    INTENT(in)  :: cldfra(KLON,KLEV), cldemi(KLON,KLEV), cldtaupd(KLON,KLEV)
    223   REAL,    INTENT(in)  :: tau_aero(KLON,KLEV,naero_grp,2)                        ! aerosol optical properties (see aeropt.F)
    224   REAL,    INTENT(in)  :: piz_aero(KLON,KLEV,naero_grp,2)                        ! aerosol optical properties (see aeropt.F)
    225   REAL,    INTENT(in)  :: cg_aero(KLON,KLEV,naero_grp,2)                         ! aerosol optical properties (see aeropt.F)
    226 !--OB
    227   REAL,    INTENT(in)  :: tau_aero_sw_rrtm(KLON,KLEV,2,NSW)                 ! aerosol optical properties RRTM
    228   REAL,    INTENT(in)  :: piz_aero_sw_rrtm(KLON,KLEV,2,NSW)                 ! aerosol optical properties RRTM
    229   REAL,    INTENT(in)  :: cg_aero_sw_rrtm(KLON,KLEV,2,NSW)                  ! aerosol optical properties RRTM
    230 ! AI
    231 !--OB fin
    232 
    233 !--C. Kleinschmitt
     91    !======================================================================
     92    ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19960719
     93    ! Objet: interface entre le modele et les rayonnements
     94    ! Arguments:
     95    !                  INPUTS
     96    ! dist----- input-R- distance astronomique terre-soleil
     97    ! rmu0----- input-R- cosinus de l'angle zenithal
     98    ! fract---- input-R- duree d'ensoleillement normalisee
     99    ! co2_ppm-- input-R- concentration du gaz carbonique (en ppm)
     100    ! paprs---- input-R- pression a inter-couche (Pa)
     101    ! pplay---- input-R- pression au milieu de couche (Pa)
     102    ! tsol----- input-R- temperature du sol (en K)
     103    ! alb1----- input-R- albedo du sol(entre 0 et 1) dans l'interval visible
     104    ! alb2----- input-R- albedo du sol(entre 0 et 1) dans l'interval proche infra-rouge   
     105    ! t-------- input-R- temperature (K)
     106    ! q-------- input-R- vapeur d'eau (en kg/kg)
     107    ! cldfra--- input-R- fraction nuageuse (entre 0 et 1)
     108    ! cldtaupd- input-R- epaisseur optique des nuages dans le visible (present-day value)
     109    ! cldemi--- input-R- emissivite des nuages dans l'IR (entre 0 et 1)
     110    ! ok_ade--- input-L- apply the Aerosol Direct Effect or not?
     111    ! ok_aie--- input-L- apply the Aerosol Indirect Effect or not?
     112    ! ok_volcan input-L- activate volcanic diags (SW heat & LW cool rate, SW & LW flux)
     113    ! flag_volc_surfstrat input-I- activate volcanic surf cooling or strato heating (or nothing)
     114    ! flag_aerosol input-I- aerosol flag from 0 to 6
     115    ! flag_aerosol_strat input-I- use stratospheric aerosols flag (0, 1, 2)
     116    ! flag_aer_feedback  input-I- activate aerosol radiative feedback (T, F)
     117    ! tau_ae, piz_ae, cg_ae input-R- aerosol optical properties (calculated in aeropt.F)
     118    ! cldtaupi  input-R- epaisseur optique des nuages dans le visible
     119    !                   calculated for pre-industrial (pi) aerosol concentrations, i.e. with smaller
     120    !                   droplet concentration, thus larger droplets, thus generally cdltaupi cldtaupd
     121    !                   it is needed for the diagnostics of the aerosol indirect radiative forcing     
     122    !
     123    !                  OUTPUTS
     124    ! heat-----output-R- echauffement atmospherique (visible) (K/jour)
     125    ! cool-----output-R- refroidissement dans l'IR (K/jour)
     126    ! albpla---output-R- albedo planetaire (entre 0 et 1)
     127    ! topsw----output-R- flux solaire net au sommet de l'atm.
     128    ! toplw----output-R- ray. IR montant au sommet de l'atmosphere
     129    ! solsw----output-R- flux solaire net a la surface
     130    ! solswfdiff----output-R- fraction de rayonnement diffus pour le flux solaire descendant a la surface
     131    ! sollw----output-R- ray. IR montant a la surface
     132    ! solswad---output-R- ray. solaire net absorbe a la surface (aerosol dir)
     133    ! topswad---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol dir)
     134    ! solswai---output-R- ray. solaire net absorbe a la surface (aerosol ind)
     135    ! topswai---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol ind)
     136    !
     137    ! heat_volc-----output-R- echauffement atmospherique  du au forcage volcanique (visible) (K/s)
     138    ! cool_volc-----output-R- refroidissement dans l'IR du au forcage volcanique (K/s)
     139    !
     140    ! ATTENTION: swai and swad have to be interpreted in the following manner:
     141    ! ---------
     142    ! ok_ade=F & ok_aie=F -both are zero
     143    ! ok_ade=T & ok_aie=F -aerosol direct forcing is F_{AD} = topsw-topswad
     144    !                        indirect is zero
     145    ! ok_ade=F & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai
     146    !                        direct is zero
     147    ! ok_ade=T & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai
     148    !                        aerosol direct forcing is F_{AD} = topswai-topswad
     149    !
     150    ! --------- RRTM: output RECMWFL
     151    ! ZEMTD (KPROMA,KLEV+1)         ; TOTAL DOWNWARD LONGWAVE EMISSIVITY
     152    ! ZEMTU (KPROMA,KLEV+1)         ; TOTAL UPWARD   LONGWAVE EMISSIVITY
     153    ! ZTRSO (KPROMA,KLEV+1)         ; TOTAL SHORTWAVE TRANSMISSIVITY
     154    ! ZTH   (KPROMA,KLEV+1)         ; HALF LEVEL TEMPERATURE
     155    ! ZCTRSO(KPROMA,2)              ; CLEAR-SKY SHORTWAVE TRANSMISSIVITY
     156    ! ZCEMTR(KPROMA,2)              ; CLEAR-SKY NET LONGWAVE EMISSIVITY
     157    ! ZTRSOD(KPROMA)                ; TOTAL-SKY SURFACE SW TRANSMISSITY
     158    ! ZLWFC (KPROMA,2)              ; CLEAR-SKY LONGWAVE FLUXES
     159    ! ZLWFT (KPROMA,KLEV+1)         ; TOTAL-SKY LONGWAVE FLUXES
     160    ! ZLWFT0(KPROMA,KLEV+1)         ; CLEAR-SKY LONGWAVE FLUXES      ! added by MPL 090109
     161    ! ZSWFC (KPROMA,2)              ; CLEAR-SKY SHORTWAVE FLUXES
     162    ! ZSWFT (KPROMA,KLEV+1)         ; TOTAL-SKY SHORTWAVE FLUXES
     163    ! ZSWFT0(KPROMA,KLEV+1)         ; CLEAR-SKY SHORTWAVE FLUXES     ! added by MPL 090109
     164    ! ZFLUX (KLON,2,KLEV+1)         ; TOTAL LW FLUXES  1=up, 2=DWN   ! added by MPL 080411
     165    ! ZFLUC (KLON,2,KLEV+1)         ; CLEAR SKY LW FLUXES            ! added by MPL 080411
     166    ! ZFSDWN(klon,KLEV+1)           ; TOTAL SW  DWN FLUXES           ! added by MPL 080411
     167    ! ZFCDWN(klon,KLEV+1)           ; CLEAR SKY SW  DWN FLUXES       ! added by MPL 080411
     168    ! ZFCCDWN(klon,KLEV+1)          ; CLEAR SKY CLEAN (NO AEROSOL) SW  DWN FLUXES      ! added by OB 211117
     169    ! ZFSUP (klon,KLEV+1)           ; TOTAL SW  UP  FLUXES           ! added by MPL 080411
     170    ! ZFCUP (klon,KLEV+1)           ; CLEAR SKY SW  UP  FLUXES       ! added by MPL 080411
     171    ! ZFCCUP (klon,KLEV+1)          ; CLEAR SKY CLEAN (NO AEROSOL) SW  UP  FLUXES      ! added by OB 211117
     172    ! ZFLCCDWN(klon,KLEV+1)         ; CLEAR SKY CLEAN (NO AEROSOL) LW  DWN FLUXES      ! added by OB 211117
     173    ! ZFLCCUP (klon,KLEV+1)         ; CLEAR SKY CLEAN (NO AEROSOL) LW  UP  FLUXES      ! added by OB 211117
     174
     175    !======================================================================
     176
     177    ! ====================================================================
     178    ! Adapte au modele de chimie INCA par Celine Deandreis & Anne Cozic -- 2009
     179    ! 1 = ZERO   
     180    ! 2 = AER total   
     181    ! 3 = NAT   
     182    ! 4 = BC   
     183    ! 5 = SO4   
     184    ! 6 = POM   
     185    ! 7 = DUST   
     186    ! 8 = SS   
     187    ! 9 = NO3   
     188    !
     189    ! ====================================================================
     190
     191    ! ==============
     192    ! DECLARATIONS
     193    ! ==============
     194    include "YOETHF.h"
     195    include "YOMCST.h"
     196    include "clesphys.h"
     197
     198    ! Input arguments
     199    REAL,    INTENT(in)  :: dist
     200    REAL,    INTENT(in)  :: rmu0(KLON), fract(KLON)
     201    REAL,    INTENT(in)  :: paprs(KLON,KLEV+1), pplay(KLON,KLEV)
     202    !albedo SB >>>
     203    ! REAL,    INTENT(in)  :: alb1(KLON), alb2(KLON), tsol(KLON)
     204    REAL,    INTENT(in)  :: tsol(KLON)
     205    REAL,    INTENT(in) :: alb_dir(KLON,NSW),alb_dif(KLON,NSW)
     206    REAL,    INTENT(in) :: SFRWL(6)
     207    !albedo SB <<<
     208    REAL,    INTENT(in)  :: t(KLON,KLEV), q(KLON,KLEV)
     209
     210    REAL, INTENT(in):: wo(:, :, :) ! dimension(KLON,KLEV, 1 or 2)
     211    ! column-density of ozone in a layer, in kilo-Dobsons
     212    ! "wo(:, :, 1)" is for the average day-night field,
     213    ! "wo(:, :, 2)" is for daylight time.
     214
     215    LOGICAL, INTENT(in)  :: ok_ade, ok_aie                                 ! switches whether to use aerosol direct (indirect) effects or not
     216    LOGICAL, INTENT(in)  :: ok_volcan                                      ! produce volcanic diags (SW/LW heat flux and rate)
     217    INTEGER, INTENT(in)  :: flag_volc_surfstrat                            ! allow to impose volcanic cooling rate at surf or heating in strato
     218    LOGICAL              :: lldebug=.false.
     219    INTEGER, INTENT(in)  :: flag_aerosol                                   ! takes value 0 (no aerosol) or 1 to 6 (aerosols)
     220    INTEGER, INTENT(in)  :: flag_aerosol_strat                             ! use stratospheric aerosols
     221    LOGICAL, INTENT(in)  :: flag_aer_feedback                              ! activate aerosol radiative feedback
     222    REAL,    INTENT(in)  :: cldfra(KLON,KLEV), cldemi(KLON,KLEV), cldtaupd(KLON,KLEV)
     223    REAL,    INTENT(in)  :: tau_aero(KLON,KLEV,naero_grp,2)                        ! aerosol optical properties (see aeropt.F)
     224    REAL,    INTENT(in)  :: piz_aero(KLON,KLEV,naero_grp,2)                        ! aerosol optical properties (see aeropt.F)
     225    REAL,    INTENT(in)  :: cg_aero(KLON,KLEV,naero_grp,2)                         ! aerosol optical properties (see aeropt.F)
     226    !--OB
     227    REAL,    INTENT(in)  :: tau_aero_sw_rrtm(KLON,KLEV,2,NSW)                 ! aerosol optical properties RRTM
     228    REAL,    INTENT(in)  :: piz_aero_sw_rrtm(KLON,KLEV,2,NSW)                 ! aerosol optical properties RRTM
     229    REAL,    INTENT(in)  :: cg_aero_sw_rrtm(KLON,KLEV,2,NSW)                  ! aerosol optical properties RRTM
     230    ! AI
     231    !--OB fin
     232
     233    !--C. Kleinschmitt
    234234#ifdef CPP_RRTM
    235   REAL,    INTENT(in)  :: tau_aero_lw_rrtm(KLON,KLEV,2,NLW)                 ! LW aerosol optical properties RRTM
     235    REAL,    INTENT(in)  :: tau_aero_lw_rrtm(KLON,KLEV,2,NLW)                 ! LW aerosol optical properties RRTM
    236236#else
    237   REAL,    INTENT(in)  :: tau_aero_lw_rrtm(KLON,KLEV,2,nbands_lw_rrtm)
     237    REAL,    INTENT(in)  :: tau_aero_lw_rrtm(KLON,KLEV,2,nbands_lw_rrtm)
    238238#endif
    239 !--C. Kleinschmitt end
    240 
    241   REAL,    INTENT(in)  :: cldtaupi(KLON,KLEV)                            ! cloud optical thickness for pre-industrial aerosol concentrations
    242   REAL,    INTENT(in)  :: qsat(klon,klev) ! Variable pour iflag_rrtm=1
    243   REAL,    INTENT(in)  :: flwc(klon,klev) ! Variable pour iflag_rrtm=1
    244   REAL,    INTENT(in)  :: fiwc(klon,klev) ! Variable pour iflag_rrtm=1
    245   REAL,    INTENT(in)  :: ref_liq(klon,klev) ! cloud droplet radius present-day from newmicro
    246   REAL,    INTENT(in)  :: ref_ice(klon,klev) ! ice crystal radius   present-day from newmicro
    247   REAL,    INTENT(in)  :: ref_liq_pi(klon,klev) ! cloud droplet radius pre-industrial from newmicro
    248   REAL,    INTENT(in)  :: ref_ice_pi(klon,klev) ! ice crystal radius   pre-industrial from newmicro
    249 
    250   CHARACTER(len=512), INTENT(in) :: namelist_ecrad_file
    251   LOGICAL, INTENT(in)  :: debut
    252 
    253 ! Output arguments
    254   REAL,    INTENT(out) :: heat(KLON,KLEV), cool(KLON,KLEV)
    255   REAL,    INTENT(out) :: heat0(KLON,KLEV), cool0(KLON,KLEV)
    256   REAL,    INTENT(out) :: heat_volc(KLON,KLEV), cool_volc(KLON,KLEV) !NL
    257   REAL,    INTENT(out) :: topsw(KLON), toplw(KLON)
    258   REAL,    INTENT(out) :: solsw(KLON), sollw(KLON), albpla(KLON), solswfdiff(KLON)
    259   REAL,    INTENT(out) :: topsw0(KLON), toplw0(KLON), solsw0(KLON), sollw0(KLON)
    260   REAL,    INTENT(out) :: sollwdown(KLON)
    261   REAL,    INTENT(out) :: swdn(KLON,kflev+1),swdn0(KLON,kflev+1), swdnc0(KLON,kflev+1)
    262   REAL,    INTENT(out) :: swup(KLON,kflev+1),swup0(KLON,kflev+1), swupc0(KLON,kflev+1)
    263   REAL,    INTENT(out) :: lwdn(KLON,kflev+1),lwdn0(KLON,kflev+1), lwdnc0(KLON,kflev+1)
    264   REAL,    INTENT(out) :: lwup(KLON,kflev+1),lwup0(KLON,kflev+1), lwupc0(KLON,kflev+1)
    265   REAL,    INTENT(out) :: topswad_aero(KLON), solswad_aero(KLON)         ! output: aerosol direct forcing at TOA and surface
    266   REAL,    INTENT(out) :: topswai_aero(KLON), solswai_aero(KLON)         ! output: aerosol indirect forcing atTOA and surface
    267   REAL,    INTENT(out) :: toplwad_aero(KLON), sollwad_aero(KLON)         ! output: LW aerosol direct forcing at TOA and surface
    268   REAL,    INTENT(out) :: toplwai_aero(KLON), sollwai_aero(KLON)         ! output: LW aerosol indirect forcing atTOA and surface
    269   REAL, DIMENSION(klon), INTENT(out)    :: topswad0_aero
    270   REAL, DIMENSION(klon), INTENT(out)    :: solswad0_aero
    271   REAL, DIMENSION(klon), INTENT(out)    :: toplwad0_aero
    272   REAL, DIMENSION(klon), INTENT(out)    :: sollwad0_aero
    273   REAL, DIMENSION(kdlon,9), INTENT(out) :: topsw_aero
    274   REAL, DIMENSION(kdlon,9), INTENT(out) :: topsw0_aero
    275   REAL, DIMENSION(kdlon,9), INTENT(out) :: solsw_aero
    276   REAL, DIMENSION(kdlon,9), INTENT(out) :: solsw0_aero
    277   REAL, DIMENSION(kdlon,3), INTENT(out) :: topswcf_aero
    278   REAL, DIMENSION(kdlon,3), INTENT(out) :: solswcf_aero
    279   REAL, DIMENSION(kdlon,kflev+1), INTENT(out) :: ZSWFT0_i
    280   REAL, DIMENSION(kdlon,kflev+1), INTENT(out) :: ZLWFT0_i
    281 
    282 ! Local variables
    283   REAL(KIND=8) ZFSUP(KDLON,KFLEV+1)
    284   REAL(KIND=8) ZFSDN(KDLON,KFLEV+1)
    285   REAL(KIND=8) ZFSUP0(KDLON,KFLEV+1)
    286   REAL(KIND=8) ZFSDN0(KDLON,KFLEV+1)
    287   REAL(KIND=8) ZFSUPC0(KDLON,KFLEV+1)
    288   REAL(KIND=8) ZFSDNC0(KDLON,KFLEV+1)
    289   REAL(KIND=8) ZFLUP(KDLON,KFLEV+1)
    290   REAL(KIND=8) ZFLDN(KDLON,KFLEV+1)
    291   REAL(KIND=8) ZFLUP0(KDLON,KFLEV+1)
    292   REAL(KIND=8) ZFLDN0(KDLON,KFLEV+1)
    293   REAL(KIND=8) ZFLUPC0(KDLON,KFLEV+1)
    294   REAL(KIND=8) ZFLDNC0(KDLON,KFLEV+1)
    295   REAL(KIND=8) zx_alpha1, zx_alpha2
    296   INTEGER k, kk, i, j, iof, nb_gr
    297   INTEGER ist,iend,ktdia,kmode
    298   REAL(KIND=8) PSCT
    299   REAL(KIND=8) PALBD(kdlon,2), PALBP(kdlon,2)
    300 !  MPL 06.01.09: pour RRTM, creation de PALBD_NEW et PALBP_NEW
    301 ! avec NSW en deuxieme dimension       
    302   REAL(KIND=8) PALBD_NEW(kdlon,NSW), PALBP_NEW(kdlon,NSW)
    303   REAL(KIND=8) PEMIS(kdlon), PDT0(kdlon), PVIEW(kdlon)
    304   REAL(KIND=8) PPSOL(kdlon), PDP(kdlon,KLEV)
    305   REAL(KIND=8) PTL(kdlon,kflev+1), PPMB(kdlon,kflev+1)
    306   REAL(KIND=8) PTAVE(kdlon,kflev)
    307   REAL(KIND=8) PWV(kdlon,kflev), PQS(kdlon,kflev)
    308 
    309   REAL(KIND=8) cloud_cover_sw(klon)
     239    !--C. Kleinschmitt end
     240
     241    REAL,    INTENT(in)  :: cldtaupi(KLON,KLEV)                            ! cloud optical thickness for pre-industrial aerosol concentrations
     242    REAL,    INTENT(in)  :: qsat(klon,klev) ! Variable pour iflag_rrtm=1
     243    REAL,    INTENT(in)  :: flwc(klon,klev) ! Variable pour iflag_rrtm=1
     244    REAL,    INTENT(in)  :: fiwc(klon,klev) ! Variable pour iflag_rrtm=1
     245    REAL,    INTENT(in)  :: ref_liq(klon,klev) ! cloud droplet radius present-day from newmicro
     246    REAL,    INTENT(in)  :: ref_ice(klon,klev) ! ice crystal radius   present-day from newmicro
     247    REAL,    INTENT(in)  :: ref_liq_pi(klon,klev) ! cloud droplet radius pre-industrial from newmicro
     248    REAL,    INTENT(in)  :: ref_ice_pi(klon,klev) ! ice crystal radius   pre-industrial from newmicro
     249
     250    CHARACTER(len=512), INTENT(in) :: namelist_ecrad_file
     251    LOGICAL, INTENT(in)  :: debut
     252
     253    ! Output arguments
     254    REAL,    INTENT(out) :: heat(KLON,KLEV), cool(KLON,KLEV)
     255    REAL,    INTENT(out) :: heat0(KLON,KLEV), cool0(KLON,KLEV)
     256    REAL,    INTENT(out) :: heat_volc(KLON,KLEV), cool_volc(KLON,KLEV) !NL
     257    REAL,    INTENT(out) :: topsw(KLON), toplw(KLON)
     258    REAL,    INTENT(out) :: solsw(KLON), sollw(KLON), albpla(KLON), solswfdiff(KLON)
     259    REAL,    INTENT(out) :: topsw0(KLON), toplw0(KLON), solsw0(KLON), sollw0(KLON)
     260    REAL,    INTENT(out) :: sollwdown(KLON)
     261    REAL,    INTENT(out) :: swdn(KLON,kflev+1),swdn0(KLON,kflev+1), swdnc0(KLON,kflev+1)
     262    REAL,    INTENT(out) :: swup(KLON,kflev+1),swup0(KLON,kflev+1), swupc0(KLON,kflev+1)
     263    REAL,    INTENT(out) :: lwdn(KLON,kflev+1),lwdn0(KLON,kflev+1), lwdnc0(KLON,kflev+1)
     264    REAL,    INTENT(out) :: lwup(KLON,kflev+1),lwup0(KLON,kflev+1), lwupc0(KLON,kflev+1)
     265    REAL,    INTENT(out) :: topswad_aero(KLON), solswad_aero(KLON)         ! output: aerosol direct forcing at TOA and surface
     266    REAL,    INTENT(out) :: topswai_aero(KLON), solswai_aero(KLON)         ! output: aerosol indirect forcing atTOA and surface
     267    REAL,    INTENT(out) :: toplwad_aero(KLON), sollwad_aero(KLON)         ! output: LW aerosol direct forcing at TOA and surface
     268    REAL,    INTENT(out) :: toplwai_aero(KLON), sollwai_aero(KLON)         ! output: LW aerosol indirect forcing atTOA and surface
     269    REAL, DIMENSION(klon), INTENT(out)    :: topswad0_aero
     270    REAL, DIMENSION(klon), INTENT(out)    :: solswad0_aero
     271    REAL, DIMENSION(klon), INTENT(out)    :: toplwad0_aero
     272    REAL, DIMENSION(klon), INTENT(out)    :: sollwad0_aero
     273    REAL, DIMENSION(kdlon,9), INTENT(out) :: topsw_aero
     274    REAL, DIMENSION(kdlon,9), INTENT(out) :: topsw0_aero
     275    REAL, DIMENSION(kdlon,9), INTENT(out) :: solsw_aero
     276    REAL, DIMENSION(kdlon,9), INTENT(out) :: solsw0_aero
     277    REAL, DIMENSION(kdlon,3), INTENT(out) :: topswcf_aero
     278    REAL, DIMENSION(kdlon,3), INTENT(out) :: solswcf_aero
     279    REAL, DIMENSION(kdlon,kflev+1), INTENT(out) :: ZSWFT0_i
     280    REAL, DIMENSION(kdlon,kflev+1), INTENT(out) :: ZLWFT0_i
     281
     282    ! Local variables
     283    REAL(KIND=8) ZFSUP(KDLON,KFLEV+1)
     284    REAL(KIND=8) ZFSDN(KDLON,KFLEV+1)
     285    REAL(KIND=8) ZFSUP0(KDLON,KFLEV+1)
     286    REAL(KIND=8) ZFSDN0(KDLON,KFLEV+1)
     287    REAL(KIND=8) ZFSUPC0(KDLON,KFLEV+1)
     288    REAL(KIND=8) ZFSDNC0(KDLON,KFLEV+1)
     289    REAL(KIND=8) ZFLUP(KDLON,KFLEV+1)
     290    REAL(KIND=8) ZFLDN(KDLON,KFLEV+1)
     291    REAL(KIND=8) ZFLUP0(KDLON,KFLEV+1)
     292    REAL(KIND=8) ZFLDN0(KDLON,KFLEV+1)
     293    REAL(KIND=8) ZFLUPC0(KDLON,KFLEV+1)
     294    REAL(KIND=8) ZFLDNC0(KDLON,KFLEV+1)
     295    REAL(KIND=8) zx_alpha1, zx_alpha2
     296    INTEGER k, kk, i, j, iof, nb_gr
     297    INTEGER ist,iend,ktdia,kmode
     298    REAL(KIND=8) PSCT
     299    REAL(KIND=8) PALBD(kdlon,2), PALBP(kdlon,2)
     300    !  MPL 06.01.09: pour RRTM, creation de PALBD_NEW et PALBP_NEW
     301    ! avec NSW en deuxieme dimension       
     302    REAL(KIND=8) PALBD_NEW(kdlon,NSW), PALBP_NEW(kdlon,NSW)
     303    REAL(KIND=8) PEMIS(kdlon), PDT0(kdlon), PVIEW(kdlon)
     304    REAL(KIND=8) PPSOL(kdlon), PDP(kdlon,KLEV)
     305    REAL(KIND=8) PTL(kdlon,kflev+1), PPMB(kdlon,kflev+1)
     306    REAL(KIND=8) PTAVE(kdlon,kflev)
     307    REAL(KIND=8) PWV(kdlon,kflev), PQS(kdlon,kflev)
     308
     309    REAL(KIND=8) cloud_cover_sw(klon)
    310310
    311311!!!!!!! Declarations specifiques pour ECRAD !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    312 ! AI 02.2021
     312    ! AI 02.2021
    313313#ifdef CPP_ECRAD
    314 ! ATTENTION les dimensions klon, kdlon ???
    315 ! INPUTS
    316   REAL, DIMENSION(kdlon,kflev+1) :: ZSWFT0_ii, ZLWFT0_ii
    317   REAL(KIND=8) ZEMISW(klon), &              ! LW emissivity inside the window region
    318                ZEMIS(klon)                  ! LW emissivity outside the window region
    319   REAL(KIND=8) ZGELAM(klon), &              ! longitudes en rad
    320                ZGEMU(klon)                  ! sin(latitude)
    321   REAL(KIND=8) ZCO2, &           ! CO2 mass mixing ratios on full levels
    322                ZCH4, &           ! CH4 mass mixing ratios on full levels
    323                ZN2O, &           ! N2O mass mixing ratios on full levels
    324                ZNO2, &           ! NO2 mass mixing ratios on full levels
    325                ZCFC11, &         ! CFC11
    326                ZCFC12, &         ! CFC12
    327                ZHCFC22, &        ! HCFC22
    328                ZCCL4, &          ! CCL4
    329                ZO2               ! O2
    330 
    331   REAL(KIND=8) ZQ_RAIN(klon,klev), &        ! Rain cloud mass mixing ratio (kg/kg) ?
    332                ZQ_SNOW(klon,klev)           ! Snow cloud mass mixing ratio (kg/kg) ?
    333   REAL(KIND=8) ZAEROSOL_OLD(KLON,6,KLEV), &  !
    334                ZAEROSOL(KLON,KLEV,naero_spc) !
    335 ! OUTPUTS
    336   REAL(KIND=8) ZFLUX_DIR(klon), &           ! Direct compt of surf flux into horizontal plane
    337                ZFLUX_DIR_CLEAR(klon), &     ! CS Direct
    338                ZFLUX_DIR_INTO_SUN(klon), &  !
    339                ZFLUX_UV(klon), &            ! UV flux
    340                ZFLUX_PAR(klon), &           ! photosynthetically active radiation similarly
    341                ZFLUX_PAR_CLEAR(klon), &     ! CS photosynthetically
    342                ZFLUX_SW_DN_TOA(klon), &     ! DN SW flux at TOA
    343                ZEMIS_OUT(klon)              ! effective broadband emissivity
    344 
    345   REAL(KIND=8) ZLWDERIVATIVE(klon,klev+1)   ! LW derivatives
    346   REAL(KIND=8) ZSWDIFFUSEBAND(klon,NSW), &  ! SW DN flux in diffuse albedo band
    347                ZSWDIRECTBAND(klon,NSW)      ! SW DN flux in direct albedo band
    348   REAL(KIND=8) SOLARIRAD
    349   REAL(KIND=8) seuilmach
    350 ! AI 10 mars 22 : Pour les tests Offline
    351   logical   :: lldebug_for_offline = .false.
    352   REAL(KIND=8) solaire_off(klon), &
    353                ZCO2_off(klon,klev), &
    354                ZCH4_off(klon,klev), &           ! CH4 mass mixing ratios on full levels
    355                ZN2O_off(klon,klev), &           ! N2O mass mixing ratios on full levels
    356                ZNO2_off(klon,klev), &           ! NO2 mass mixing ratios on full levels
    357                ZCFC11_off(klon,klev), &         ! CFC11
    358                ZCFC12_off(klon,klev), &         ! CFC12
    359                ZHCFC22_off(klon,klev), &        ! HCFC22
    360                ZCCL4_off(klon,klev), &          ! CCL4
    361                ZO2_off(klon,klev)               ! O2#endif
     314    ! ATTENTION les dimensions klon, kdlon ???
     315    ! INPUTS
     316    REAL, DIMENSION(kdlon,kflev+1) :: ZSWFT0_ii, ZLWFT0_ii
     317    REAL(KIND=8) ZEMISW(klon), &              ! LW emissivity inside the window region
     318         ZEMIS(klon)                  ! LW emissivity outside the window region
     319    REAL(KIND=8) ZGELAM(klon), &              ! longitudes en rad
     320         ZGEMU(klon)                  ! sin(latitude)
     321    REAL(KIND=8) ZCO2, &           ! CO2 mass mixing ratios on full levels
     322         ZCH4, &           ! CH4 mass mixing ratios on full levels
     323         ZN2O, &           ! N2O mass mixing ratios on full levels
     324         ZNO2, &           ! NO2 mass mixing ratios on full levels
     325         ZCFC11, &         ! CFC11
     326         ZCFC12, &         ! CFC12
     327         ZHCFC22, &        ! HCFC22
     328         ZCCL4, &          ! CCL4
     329         ZO2               ! O2
     330
     331    REAL(KIND=8) ZQ_RAIN(klon,klev), &        ! Rain cloud mass mixing ratio (kg/kg) ?
     332         ZQ_SNOW(klon,klev)           ! Snow cloud mass mixing ratio (kg/kg) ?
     333    REAL(KIND=8) ZAEROSOL_OLD(KLON,6,KLEV), &  !
     334         ZAEROSOL(KLON,KLEV,naero_spc) !
     335    ! OUTPUTS
     336    REAL(KIND=8) ZFLUX_DIR(klon), &           ! Direct compt of surf flux into horizontal plane
     337         ZFLUX_DIR_CLEAR(klon), &     ! CS Direct
     338         ZFLUX_DIR_INTO_SUN(klon), &  !
     339         ZFLUX_UV(klon), &            ! UV flux
     340         ZFLUX_PAR(klon), &           ! photosynthetically active radiation similarly
     341         ZFLUX_PAR_CLEAR(klon), &     ! CS photosynthetically
     342         ZFLUX_SW_DN_TOA(klon), &     ! DN SW flux at TOA
     343         ZEMIS_OUT(klon)              ! effective broadband emissivity
     344
     345    REAL(KIND=8) ZLWDERIVATIVE(klon,klev+1)   ! LW derivatives
     346    REAL(KIND=8) ZSWDIFFUSEBAND(klon,NSW), &  ! SW DN flux in diffuse albedo band
     347         ZSWDIRECTBAND(klon,NSW)      ! SW DN flux in direct albedo band
     348    REAL(KIND=8) SOLARIRAD
     349    REAL(KIND=8) seuilmach
     350    ! AI 10 mars 22 : Pour les tests Offline
     351    logical   :: lldebug_for_offline = .false.
     352    REAL(KIND=8) solaire_off(klon), &
     353         ZCO2_off(klon,klev), &
     354         ZCH4_off(klon,klev), &           ! CH4 mass mixing ratios on full levels
     355         ZN2O_off(klon,klev), &           ! N2O mass mixing ratios on full levels
     356         ZNO2_off(klon,klev), &           ! NO2 mass mixing ratios on full levels
     357         ZCFC11_off(klon,klev), &         ! CFC11
     358         ZCFC12_off(klon,klev), &         ! CFC12
     359         ZHCFC22_off(klon,klev), &        ! HCFC22
     360         ZCCL4_off(klon,klev), &          ! CCL4
     361         ZO2_off(klon,klev)               ! O2#endif
    362362#endif
    363363!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    364364
    365   REAL(kind=8) POZON(kdlon, kflev, size(wo, 3)) ! mass fraction of ozone
    366   ! "POZON(:, :, 1)" is for the average day-night field,
    367   ! "POZON(:, :, 2)" is for daylight time.
     365    REAL(kind=8) POZON(kdlon, kflev, size(wo, 3)) ! mass fraction of ozone
     366    ! "POZON(:, :, 1)" is for the average day-night field,
     367    ! "POZON(:, :, 2)" is for daylight time.
    368368!!!!! Modif MPL 6.01.09 avec RRTM, on passe de 5 a 6 
    369   REAL(KIND=8) PAER(kdlon,kflev,6)
    370   REAL(KIND=8) PCLDLD(kdlon,kflev)
    371   REAL(KIND=8) PCLDLU(kdlon,kflev)
    372   REAL(KIND=8) PCLDSW(kdlon,kflev)
    373   REAL(KIND=8) PTAU(kdlon,2,kflev)
    374   REAL(KIND=8) POMEGA(kdlon,2,kflev)
    375   REAL(KIND=8) PCG(kdlon,2,kflev)
    376   REAL(KIND=8) zfract(kdlon), zrmu0(kdlon), zdist
    377   REAL(KIND=8) zheat(kdlon,kflev), zcool(kdlon,kflev)
    378   REAL(KIND=8) zheat0(kdlon,kflev), zcool0(kdlon,kflev)
    379   REAL(KIND=8) zheat_volc(kdlon,kflev), zcool_volc(kdlon,kflev) !NL
    380   REAL(KIND=8) ztopsw(kdlon), ztoplw(kdlon)
    381   REAL(KIND=8) zsolsw(kdlon), zsollw(kdlon), zalbpla(kdlon), zsolswfdiff(kdlon)
    382   REAL(KIND=8) zsollwdown(kdlon)
    383   REAL(KIND=8) ztopsw0(kdlon), ztoplw0(kdlon)
    384   REAL(KIND=8) zsolsw0(kdlon), zsollw0(kdlon)
    385   REAL(KIND=8) zznormcp
    386   REAL(KIND=8) tauaero(kdlon,kflev,naero_grp,2)                     ! aer opt properties
    387   REAL(KIND=8) pizaero(kdlon,kflev,naero_grp,2)
    388   REAL(KIND=8) cgaero(kdlon,kflev,naero_grp,2)
    389   REAL(KIND=8) PTAUA(kdlon,2,kflev)                         ! present-day value of cloud opt thickness (PTAU is pre-industrial value), local use
    390   REAL(KIND=8) POMEGAA(kdlon,2,kflev)                       ! dito for single scatt albedo
    391   REAL(KIND=8) ztopswadaero(kdlon), zsolswadaero(kdlon)     ! Aerosol direct forcing at TOAand surface
    392   REAL(KIND=8) ztopswad0aero(kdlon), zsolswad0aero(kdlon)   ! Aerosol direct forcing at TOAand surface
    393   REAL(KIND=8) ztopswaiaero(kdlon), zsolswaiaero(kdlon)     ! dito, indirect
    394 !--NL
    395   REAL(KIND=8) zswadaero(kdlon,kflev+1)                     ! SW Aerosol direct forcing
    396   REAL(KIND=8) zlwadaero(kdlon,kflev+1)                     ! LW Aerosol direct forcing
    397   REAL(KIND=8) volmip_solsw(kdlon)                          ! SW clear sky in the case of VOLMIP
    398 !-LW by CK
    399   REAL(KIND=8) ztoplwadaero(kdlon), zsollwadaero(kdlon)     ! LW Aerosol direct forcing at TOAand surface
    400   REAL(KIND=8) ztoplwad0aero(kdlon), zsollwad0aero(kdlon)   ! LW Aerosol direct forcing at TOAand surface
    401   REAL(KIND=8) ztoplwaiaero(kdlon), zsollwaiaero(kdlon)     ! dito, indirect
    402 !-end
    403   REAL(KIND=8) ztopsw_aero(kdlon,9), ztopsw0_aero(kdlon,9)
    404   REAL(KIND=8) zsolsw_aero(kdlon,9), zsolsw0_aero(kdlon,9)
    405   REAL(KIND=8) ztopswcf_aero(kdlon,3), zsolswcf_aero(kdlon,3)     
    406 ! real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2 deje declare dans physiq.F MPL 20130618
    407 !MPL input supplementaires pour RECMWFL
    408 ! flwc, fiwc = Liquid Water Content & Ice Water Content (kg/kg)
    409   REAL(KIND=8) GEMU(klon)
    410 !MPL input RECMWFL:
    411 ! Tableaux aux niveaux inverses pour respecter convention Arpege
    412   REAL(KIND=8) ref_liq_i(klon,klev) ! cloud droplet radius present-day from newmicro (inverted)
    413   REAL(KIND=8) ref_ice_i(klon,klev) ! ice crystal radius present-day from newmicro (inverted)
    414 !--OB
    415   REAL(KIND=8) ref_liq_pi_i(klon,klev) ! cloud droplet radius pre-industrial from newmicro (inverted)
    416   REAL(KIND=8) ref_ice_pi_i(klon,klev) ! ice crystal radius pre-industrial from newmicro (inverted)
    417 !--end OB
    418   REAL(KIND=8) paprs_i(klon,klev+1)
    419   REAL(KIND=8) pplay_i(klon,klev)
    420   REAL(KIND=8) cldfra_i(klon,klev)
    421   REAL(KIND=8) POZON_i(kdlon,kflev, size(wo, 3)) ! mass fraction of ozone
    422   ! "POZON(:, :, 1)" is for the average day-night field,
    423   ! "POZON(:, :, 2)" is for daylight time.
     369    REAL(KIND=8) PAER(kdlon,kflev,6)
     370    REAL(KIND=8) PCLDLD(kdlon,kflev)
     371    REAL(KIND=8) PCLDLU(kdlon,kflev)
     372    REAL(KIND=8) PCLDSW(kdlon,kflev)
     373    REAL(KIND=8) PTAU(kdlon,2,kflev)
     374    REAL(KIND=8) POMEGA(kdlon,2,kflev)
     375    REAL(KIND=8) PCG(kdlon,2,kflev)
     376    REAL(KIND=8) zfract(kdlon), zrmu0(kdlon), zdist
     377    REAL(KIND=8) zheat(kdlon,kflev), zcool(kdlon,kflev)
     378    REAL(KIND=8) zheat0(kdlon,kflev), zcool0(kdlon,kflev)
     379    REAL(KIND=8) zheat_volc(kdlon,kflev), zcool_volc(kdlon,kflev) !NL
     380    REAL(KIND=8) ztopsw(kdlon), ztoplw(kdlon)
     381    REAL(KIND=8) zsolsw(kdlon), zsollw(kdlon), zalbpla(kdlon), zsolswfdiff(kdlon)
     382    REAL(KIND=8) zsollwdown(kdlon)
     383    REAL(KIND=8) ztopsw0(kdlon), ztoplw0(kdlon)
     384    REAL(KIND=8) zsolsw0(kdlon), zsollw0(kdlon)
     385    REAL(KIND=8) zznormcp
     386    REAL(KIND=8) tauaero(kdlon,kflev,naero_grp,2)                     ! aer opt properties
     387    REAL(KIND=8) pizaero(kdlon,kflev,naero_grp,2)
     388    REAL(KIND=8) cgaero(kdlon,kflev,naero_grp,2)
     389    REAL(KIND=8) PTAUA(kdlon,2,kflev)                         ! present-day value of cloud opt thickness (PTAU is pre-industrial value), local use
     390    REAL(KIND=8) POMEGAA(kdlon,2,kflev)                       ! dito for single scatt albedo
     391    REAL(KIND=8) ztopswadaero(kdlon), zsolswadaero(kdlon)     ! Aerosol direct forcing at TOAand surface
     392    REAL(KIND=8) ztopswad0aero(kdlon), zsolswad0aero(kdlon)   ! Aerosol direct forcing at TOAand surface
     393    REAL(KIND=8) ztopswaiaero(kdlon), zsolswaiaero(kdlon)     ! dito, indirect
     394    !--NL
     395    REAL(KIND=8) zswadaero(kdlon,kflev+1)                     ! SW Aerosol direct forcing
     396    REAL(KIND=8) zlwadaero(kdlon,kflev+1)                     ! LW Aerosol direct forcing
     397    REAL(KIND=8) volmip_solsw(kdlon)                          ! SW clear sky in the case of VOLMIP
     398    !-LW by CK
     399    REAL(KIND=8) ztoplwadaero(kdlon), zsollwadaero(kdlon)     ! LW Aerosol direct forcing at TOAand surface
     400    REAL(KIND=8) ztoplwad0aero(kdlon), zsollwad0aero(kdlon)   ! LW Aerosol direct forcing at TOAand surface
     401    REAL(KIND=8) ztoplwaiaero(kdlon), zsollwaiaero(kdlon)     ! dito, indirect
     402    !-end
     403    REAL(KIND=8) ztopsw_aero(kdlon,9), ztopsw0_aero(kdlon,9)
     404    REAL(KIND=8) zsolsw_aero(kdlon,9), zsolsw0_aero(kdlon,9)
     405    REAL(KIND=8) ztopswcf_aero(kdlon,3), zsolswcf_aero(kdlon,3)     
     406    ! real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2 deje declare dans physiq.F MPL 20130618
     407    !MPL input supplementaires pour RECMWFL
     408    ! flwc, fiwc = Liquid Water Content & Ice Water Content (kg/kg)
     409    REAL(KIND=8) GEMU(klon)
     410    !MPL input RECMWFL:
     411    ! Tableaux aux niveaux inverses pour respecter convention Arpege
     412    REAL(KIND=8) ref_liq_i(klon,klev) ! cloud droplet radius present-day from newmicro (inverted)
     413    REAL(KIND=8) ref_ice_i(klon,klev) ! ice crystal radius present-day from newmicro (inverted)
     414    !--OB
     415    REAL(KIND=8) ref_liq_pi_i(klon,klev) ! cloud droplet radius pre-industrial from newmicro (inverted)
     416    REAL(KIND=8) ref_ice_pi_i(klon,klev) ! ice crystal radius pre-industrial from newmicro (inverted)
     417    !--end OB
     418    REAL(KIND=8) paprs_i(klon,klev+1)
     419    REAL(KIND=8) pplay_i(klon,klev)
     420    REAL(KIND=8) cldfra_i(klon,klev)
     421    REAL(KIND=8) POZON_i(kdlon,kflev, size(wo, 3)) ! mass fraction of ozone
     422    ! "POZON(:, :, 1)" is for the average day-night field,
     423    ! "POZON(:, :, 2)" is for daylight time.
    424424!!!!! Modif MPL 6.01.09 avec RRTM, on passe de 5 a 6     
    425   REAL(KIND=8) PAER_i(kdlon,kflev,6)
    426   REAL(KIND=8) PDP_i(klon,klev)
    427   REAL(KIND=8) t_i(klon,klev),q_i(klon,klev),qsat_i(klon,klev)
    428   REAL(KIND=8) flwc_i(klon,klev),fiwc_i(klon,klev)
    429 !MPL output RECMWFL:
    430   REAL(KIND=8) ZEMTD (klon,klev+1),ZEMTD_i (klon,klev+1)       
    431   REAL(KIND=8) ZEMTU (klon,klev+1),ZEMTU_i (klon,klev+1)     
    432   REAL(KIND=8) ZTRSO (klon,klev+1),ZTRSO_i (klon,klev+1)   
    433   REAL(KIND=8) ZTH   (klon,klev+1),ZTH_i   (klon,klev+1)   
    434   REAL(KIND=8) ZCTRSO(klon,2)       
    435   REAL(KIND=8) ZCEMTR(klon,2)     
    436   REAL(KIND=8) ZTRSOD(klon)       
    437   REAL(KIND=8) ZLWFC (klon,2)     
    438   REAL(KIND=8) ZLWFT (klon,klev+1),ZLWFT_i (klon,klev+1)   
    439   REAL(KIND=8) ZSWFC (klon,2)     
    440   REAL(KIND=8) ZSWFT (klon,klev+1),ZSWFT_i (klon,klev+1)
    441   REAL(KIND=8) ZFLUCDWN_i(klon,klev+1),ZFLUCUP_i(klon,klev+1)
    442   REAL(KIND=8) PPIZA_TOT(klon,klev,NSW)
    443   REAL(KIND=8) PCGA_TOT(klon,klev,NSW)
    444   REAL(KIND=8) PTAU_TOT(klon,klev,NSW)
    445   REAL(KIND=8) PPIZA_NAT(klon,klev,NSW)
    446   REAL(KIND=8) PCGA_NAT(klon,klev,NSW)
    447   REAL(KIND=8) PTAU_NAT(klon,klev,NSW)
     425    REAL(KIND=8) PAER_i(kdlon,kflev,6)
     426    REAL(KIND=8) PDP_i(klon,klev)
     427    REAL(KIND=8) t_i(klon,klev),q_i(klon,klev),qsat_i(klon,klev)
     428    REAL(KIND=8) flwc_i(klon,klev),fiwc_i(klon,klev)
     429    !MPL output RECMWFL:
     430    REAL(KIND=8) ZEMTD (klon,klev+1),ZEMTD_i (klon,klev+1)       
     431    REAL(KIND=8) ZEMTU (klon,klev+1),ZEMTU_i (klon,klev+1)     
     432    REAL(KIND=8) ZTRSO (klon,klev+1),ZTRSO_i (klon,klev+1)   
     433    REAL(KIND=8) ZTH   (klon,klev+1),ZTH_i   (klon,klev+1)   
     434    REAL(KIND=8) ZCTRSO(klon,2)       
     435    REAL(KIND=8) ZCEMTR(klon,2)     
     436    REAL(KIND=8) ZTRSOD(klon)       
     437    REAL(KIND=8) ZLWFC (klon,2)     
     438    REAL(KIND=8) ZLWFT (klon,klev+1),ZLWFT_i (klon,klev+1)   
     439    REAL(KIND=8) ZSWFC (klon,2)     
     440    REAL(KIND=8) ZSWFT (klon,klev+1),ZSWFT_i (klon,klev+1)
     441    REAL(KIND=8) ZFLUCDWN_i(klon,klev+1),ZFLUCUP_i(klon,klev+1)
     442    REAL(KIND=8) PPIZA_TOT(klon,klev,NSW)
     443    REAL(KIND=8) PCGA_TOT(klon,klev,NSW)
     444    REAL(KIND=8) PTAU_TOT(klon,klev,NSW)
     445    REAL(KIND=8) PPIZA_NAT(klon,klev,NSW)
     446    REAL(KIND=8) PCGA_NAT(klon,klev,NSW)
     447    REAL(KIND=8) PTAU_NAT(klon,klev,NSW)
    448448#ifdef CPP_RRTM
    449   REAL(KIND=8) PTAU_LW_TOT(klon,klev,NLW)
    450   REAL(KIND=8) PTAU_LW_NAT(klon,klev,NLW)
     449    REAL(KIND=8) PTAU_LW_TOT(klon,klev,NLW)
     450    REAL(KIND=8) PTAU_LW_NAT(klon,klev,NLW)
    451451#endif
    452   REAL(KIND=8) PSFSWDIR(klon,NSW)
    453   REAL(KIND=8) PSFSWDIF(klon,NSW)
    454   REAL(KIND=8) PFSDNN(klon)
    455   REAL(KIND=8) PFSDNV(klon)
    456 !MPL On ne redefinit pas les tableaux ZFLUX,ZFLUC,
    457 !MPL ZFSDWN,ZFCDWN,ZFSUP,ZFCUP car ils existent deja
    458 !MPL sous les noms de ZFLDN,ZFLDN0,ZFLUP,ZFLUP0,
    459 !MPL ZFSDN,ZFSDN0,ZFSUP,ZFSUP0
    460   REAL(KIND=8) ZFLUX_i (klon,2,klev+1)
    461   REAL(KIND=8) ZFLUC_i (klon,2,klev+1)
    462   REAL(KIND=8) ZFSDWN_i (klon,klev+1)
    463   REAL(KIND=8) ZFCDWN_i (klon,klev+1)
    464   REAL(KIND=8) ZFCCDWN_i (klon,klev+1)
    465   REAL(KIND=8) ZFSUP_i (klon,klev+1)
    466   REAL(KIND=8) ZFCUP_i (klon,klev+1)
    467   REAL(KIND=8) ZFCCUP_i (klon,klev+1)
    468   REAL(KIND=8) ZFLCCDWN_i (klon,klev+1)
    469   REAL(KIND=8) ZFLCCUP_i (klon,klev+1)
    470 ! 3 lignes suivantes a activer pour CCMVAL (MPL 20100412)
    471 !      REAL(KIND=8) RSUN(3,2)
    472 !      REAL(KIND=8) SUN(3)
    473 !      REAL(KIND=8) SUN_FRACT(2)
    474   REAL, PARAMETER:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
    475   CHARACTER (LEN=80) :: abort_message
    476   CHARACTER (LEN=80) :: modname='radlwsw_m'
    477 
    478   REAL zdir, zdif
    479 
    480 ! =========  INITIALISATIONS ==============================================
    481  IF (lldebug) THEN
    482   print*,'Entree dans radlwsw '
    483   print*,'************* INITIALISATIONS *****************************'
    484   print*,'klon, kdlon, klev, kflev =',klon, kdlon, klev, kflev
    485  ENDIF
    486 
    487   CALL assert(size(wo, 1) == klon, size(wo, 2) == klev, "radlwsw wo")
    488  
    489   ist=1
    490   iend=klon
    491   ktdia=1
    492   kmode=ist
    493 ! Aeros
    494   tauaero(:,:,:,:)=0.
    495   pizaero(:,:,:,:)=0.
    496   cgaero(:,:,:,:)=0.
    497 !  lldebug=.FALSE.
    498 
    499   ztopsw_aero(:,:)  = 0. !ym missing init : warning : not initialized in SW_AEROAR4
    500   ztopsw0_aero(:,:) = 0. !ym missing init : warning : not initialized in SW_AEROAR4
    501   zsolsw_aero(:,:)  = 0. !ym missing init : warning : not initialized in SW_AEROAR4
    502   zsolsw0_aero(:,:) = 0. !ym missing init : warning : not initialized in SW_AEROAR4
    503 
    504   ZTOPSWADAERO(:)  = 0. !ym missing init
    505   ZSOLSWADAERO(:)  = 0. !ym missing init
    506   ZTOPSWAD0AERO(:) = 0. !ym missing init
    507   ZSOLSWAD0AERO(:) = 0. !ym missing init
    508   ZTOPSWAIAERO(:)  = 0. !ym missing init
    509   ZSOLSWAIAERO(:)  = 0. !ym missing init 
    510   ZTOPSWCF_AERO(:,:)= 0.!ym missing init 
    511   ZSOLSWCF_AERO(:,:) =0. !ym missing init 
    512 
    513   !
    514 ! AI 02.2021
     452    REAL(KIND=8) PSFSWDIR(klon,NSW)
     453    REAL(KIND=8) PSFSWDIF(klon,NSW)
     454    REAL(KIND=8) PFSDNN(klon)
     455    REAL(KIND=8) PFSDNV(klon)
     456    !MPL On ne redefinit pas les tableaux ZFLUX,ZFLUC,
     457    !MPL ZFSDWN,ZFCDWN,ZFSUP,ZFCUP car ils existent deja
     458    !MPL sous les noms de ZFLDN,ZFLDN0,ZFLUP,ZFLUP0,
     459    !MPL ZFSDN,ZFSDN0,ZFSUP,ZFSUP0
     460    REAL(KIND=8) ZFLUX_i (klon,2,klev+1)
     461    REAL(KIND=8) ZFLUC_i (klon,2,klev+1)
     462    REAL(KIND=8) ZFSDWN_i (klon,klev+1)
     463    REAL(KIND=8) ZFCDWN_i (klon,klev+1)
     464    REAL(KIND=8) ZFCCDWN_i (klon,klev+1)
     465    REAL(KIND=8) ZFSUP_i (klon,klev+1)
     466    REAL(KIND=8) ZFCUP_i (klon,klev+1)
     467    REAL(KIND=8) ZFCCUP_i (klon,klev+1)
     468    REAL(KIND=8) ZFLCCDWN_i (klon,klev+1)
     469    REAL(KIND=8) ZFLCCUP_i (klon,klev+1)
     470    ! 3 lignes suivantes a activer pour CCMVAL (MPL 20100412)
     471    !      REAL(KIND=8) RSUN(3,2)
     472    !      REAL(KIND=8) SUN(3)
     473    !      REAL(KIND=8) SUN_FRACT(2)
     474    REAL, PARAMETER:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
     475    CHARACTER (LEN=80) :: abort_message
     476    CHARACTER (LEN=80) :: modname='radlwsw_m'
     477
     478    REAL zdir, zdif
     479
     480    ! =========  INITIALISATIONS ==============================================
     481    IF (lldebug) THEN
     482       print*,'Entree dans radlwsw '
     483       print*,'************* INITIALISATIONS *****************************'
     484       print*,'klon, kdlon, klev, kflev =',klon, kdlon, klev, kflev
     485    ENDIF
     486
     487    CALL assert(size(wo, 1) == klon, size(wo, 2) == klev, "radlwsw wo")
     488
     489    ist=1
     490    iend=klon
     491    ktdia=1
     492    kmode=ist
     493    ! Aeros
     494    tauaero(:,:,:,:)=0.
     495    pizaero(:,:,:,:)=0.
     496    cgaero(:,:,:,:)=0.
     497    !  lldebug=.FALSE.
     498
     499    ztopsw_aero(:,:)  = 0. !ym missing init : warning : not initialized in SW_AEROAR4
     500    ztopsw0_aero(:,:) = 0. !ym missing init : warning : not initialized in SW_AEROAR4
     501    zsolsw_aero(:,:)  = 0. !ym missing init : warning : not initialized in SW_AEROAR4
     502    zsolsw0_aero(:,:) = 0. !ym missing init : warning : not initialized in SW_AEROAR4
     503
     504    ZTOPSWADAERO(:)  = 0. !ym missing init
     505    ZSOLSWADAERO(:)  = 0. !ym missing init
     506    ZTOPSWAD0AERO(:) = 0. !ym missing init
     507    ZSOLSWAD0AERO(:) = 0. !ym missing init
     508    ZTOPSWAIAERO(:)  = 0. !ym missing init
     509    ZSOLSWAIAERO(:)  = 0. !ym missing init 
     510    ZTOPSWCF_AERO(:,:)= 0.!ym missing init 
     511    ZSOLSWCF_AERO(:,:) =0. !ym missing init 
     512
     513    !
     514    ! AI 02.2021
    515515#ifdef CPP_ECRAD
    516   ZEMIS = 1.0
    517   ZEMISW = 1.0
    518   ZGELAM = longitude
    519   ZGEMU = sin(latitude)
    520   ZCO2 = RCO2
    521   ZCH4 = RCH4
    522   ZN2O = RN2O
    523   ZNO2 = 0.0
    524   ZCFC11 = RCFC11
    525   ZCFC12 = RCFC12
    526   ZHCFC22 = 0.0
    527   ZO2 = 0.0
    528   ZCCL4 = 0.0
    529   ZQ_RAIN = 0.0
    530   ZQ_SNOW = 0.0
    531   ZAEROSOL_OLD = 0.0
    532   ZAEROSOL = 0.0
    533   seuilmach=tiny(seuilmach)
     516    ZEMIS = 1.0
     517    ZEMISW = 1.0
     518    ZGELAM = longitude
     519    ZGEMU = sin(latitude)
     520    ZCO2 = RCO2
     521    ZCH4 = RCH4
     522    ZN2O = RN2O
     523    ZNO2 = 0.0
     524    ZCFC11 = RCFC11
     525    ZCFC12 = RCFC12
     526    ZHCFC22 = 0.0
     527    ZO2 = 0.0
     528    ZCCL4 = 0.0
     529    ZQ_RAIN = 0.0
     530    ZQ_SNOW = 0.0
     531    ZAEROSOL_OLD = 0.0
     532    ZAEROSOL = 0.0
     533    seuilmach=tiny(seuilmach)
    534534#endif
    535535
    536   !-------------------------------------------
    537   nb_gr = KLON / kdlon
    538   IF (nb_gr*kdlon .NE. KLON) THEN
    539       PRINT*, "kdlon mauvais:", KLON, kdlon, nb_gr
    540       call abort_physic("radlwsw", "", 1)
    541   ENDIF
    542   IF (kflev .NE. KLEV) THEN
    543       PRINT*, "kflev differe de KLEV, kflev, KLEV"
    544       call abort_physic("radlwsw", "", 1)
    545   ENDIF
    546   !-------------------------------------------
    547   DO k = 1, KLEV
    548     DO i = 1, KLON
    549       heat(i,k)=0.
    550       cool(i,k)=0.
    551       heat_volc(i,k)=0. !NL
    552       cool_volc(i,k)=0. !NL
    553       heat0(i,k)=0.
    554       cool0(i,k)=0.
     536    !-------------------------------------------
     537    nb_gr = KLON / kdlon
     538    IF (nb_gr*kdlon .NE. KLON) THEN
     539       PRINT*, "kdlon mauvais:", KLON, kdlon, nb_gr
     540       call abort_physic("radlwsw", "", 1)
     541    ENDIF
     542    IF (kflev .NE. KLEV) THEN
     543       PRINT*, "kflev differe de KLEV, kflev, KLEV"
     544       call abort_physic("radlwsw", "", 1)
     545    ENDIF
     546    !-------------------------------------------
     547    DO k = 1, KLEV
     548       DO i = 1, KLON
     549          heat(i,k)=0.
     550          cool(i,k)=0.
     551          heat_volc(i,k)=0. !NL
     552          cool_volc(i,k)=0. !NL
     553          heat0(i,k)=0.
     554          cool0(i,k)=0.
     555       ENDDO
    555556    ENDDO
    556   ENDDO
    557   !
    558   zdist = dist
    559   !
    560   PSCT = solaire/zdist/zdist
    561 
    562   IF (type_trac == 'repr') THEN
    563 #ifdef REPROBUS
    564     IF (iflag_rrtm==0) THEN
    565       IF (ok_SUNTIME) PSCT = solaireTIME/zdist/zdist
    566       print*,'Constante solaire: ',PSCT*zdist*zdist
    567     ENDIF
    568 #endif
    569   ENDIF
    570 
    571  IF (lldebug) THEN
    572   print*,'************** Debut boucle de 1 a ', nb_gr
    573  ENDIF
    574 
    575   DO j = 1, nb_gr
    576     iof = kdlon*(j-1)
    577     DO i = 1, kdlon
    578       zfract(i) = fract(iof+i)
    579       zrmu0(i) = rmu0(iof+i)
    580 
    581 
    582       IF (iflag_rrtm==0) THEN
    583 !     Albedo
    584         PALBD(i,1)=alb_dif(iof+i,1)
    585         PALBD(i,2)=alb_dif(iof+i,2)
    586         PALBP(i,1)=alb_dir(iof+i,1)
    587         PALBP(i,2)=alb_dir(iof+i,2)
    588 ! AI 02.2021 cas iflag_rrtm=1 et 2
    589        ELSEIF (iflag_rrtm==1.OR.iflag_rrtm==2) THEN
    590         DO kk=1,NSW
    591           PALBD_NEW(i,kk)=alb_dif(iof+i,kk)
    592           PALBP_NEW(i,kk)=alb_dir(iof+i,kk)
    593         ENDDO
    594 !
    595       ENDIF
    596 !albedo SB <<<
    597 
    598       PEMIS(i) = 1.0    !!!!! A REVOIR (MPL)
    599       PVIEW(i) = 1.66
    600       PPSOL(i) = paprs(iof+i,1)
    601       zx_alpha1 = (paprs(iof+i,1)-pplay(iof+i,2))/(pplay(iof+i,1)-pplay(iof+i,2))
    602       zx_alpha2 = 1.0 - zx_alpha1
    603       PTL(i,1) = t(iof+i,1) * zx_alpha1 + t(iof+i,2) * zx_alpha2
    604       PTL(i,KLEV+1) = t(iof+i,KLEV)
    605       PDT0(i) = tsol(iof+i) - PTL(i,1)
    606     ENDDO
    607     DO k = 2, kflev
    608       DO i = 1, kdlon
    609         PTL(i,k) = (t(iof+i,k)+t(iof+i,k-1))*0.5
    610       ENDDO
    611     ENDDO
    612     DO k = 1, kflev
    613       DO i = 1, kdlon
    614         PDP(i,k) = paprs(iof+i,k)-paprs(iof+i,k+1)
    615         PTAVE(i,k) = t(iof+i,k)
    616         PWV(i,k) = MAX (q(iof+i,k), 1.0e-12)
    617         PQS(i,k) = PWV(i,k)
    618 !       Confert from  column density of ozone in a cell, in kDU, to a mass fraction
    619         POZON(i,k, :) = wo(iof+i, k, :) * RG * dobson_u * 1e3 &
    620              / (paprs(iof+i, k) - paprs(iof+i, k+1))
    621 !       A activer pour CCMVAL on prend l'ozone impose (MPL 07042010)
    622 !       POZON(i,k,:) = wo(i,k,:) 
    623 !       print *,'RADLWSW: POZON',k, POZON(i,k,1)
    624         PCLDLD(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)
    625         PCLDLU(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)
    626         PCLDSW(i,k) = cldfra(iof+i,k)
    627         PTAU(i,1,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! 1e-12 serait instable
    628         PTAU(i,2,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! pour 32-bit machines
    629         POMEGA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAU(i,1,k))
    630         POMEGA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAU(i,2,k))
    631         PCG(i,1,k) = 0.865
    632         PCG(i,2,k) = 0.910
    633         !-
    634         ! Introduced for aerosol indirect forcings.
    635         ! The following values use the cloud optical thickness calculated from
    636         ! present-day aerosol concentrations whereas the quantities without the
    637         ! "A" at the end are for pre-industial (natural-only) aerosol concentrations
    638         !
    639         PTAUA(i,1,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! 1e-12 serait instable
    640         PTAUA(i,2,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! pour 32-bit machines
    641         POMEGAA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAUA(i,1,k))
    642         POMEGAA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAUA(i,2,k))
    643       ENDDO
    644     ENDDO
     557    !
     558    zdist = dist
     559    !
     560    PSCT = solaire/zdist/zdist
    645561
    646562    IF (type_trac == 'repr') THEN
    647563#ifdef REPROBUS
    648        ndimozon = size(wo, 3)
    649        CALL RAD_INTERACTIF(POZON,iof)
     564       IF (iflag_rrtm==0) THEN
     565          IF (ok_SUNTIME) PSCT = solaireTIME/zdist/zdist
     566          print*,'Constante solaire: ',PSCT*zdist*zdist
     567       ENDIF
    650568#endif
    651569    ENDIF
    652     !
    653     DO k = 1, kflev+1
    654       DO i = 1, kdlon
    655         PPMB(i,k) = paprs(iof+i,k)/100.0
    656       ENDDO
    657     ENDDO
    658     !
     570
     571    IF (lldebug) THEN
     572       print*,'************** Debut boucle de 1 a ', nb_gr
     573    ENDIF
     574
     575    DO j = 1, nb_gr
     576       iof = kdlon*(j-1)
     577       DO i = 1, kdlon
     578          zfract(i) = fract(iof+i)
     579          zrmu0(i) = rmu0(iof+i)
     580
     581
     582          IF (iflag_rrtm==0) THEN
     583             !     Albedo
     584             PALBD(i,1)=alb_dif(iof+i,1)
     585             PALBD(i,2)=alb_dif(iof+i,2)
     586             PALBP(i,1)=alb_dir(iof+i,1)
     587             PALBP(i,2)=alb_dir(iof+i,2)
     588             ! AI 02.2021 cas iflag_rrtm=1 et 2
     589          ELSEIF (iflag_rrtm==1.OR.iflag_rrtm==2) THEN
     590             DO kk=1,NSW
     591                PALBD_NEW(i,kk)=alb_dif(iof+i,kk)
     592                PALBP_NEW(i,kk)=alb_dir(iof+i,kk)
     593             ENDDO
     594             !
     595          ENDIF
     596          !albedo SB <<<
     597
     598          PEMIS(i) = 1.0    !!!!! A REVOIR (MPL)
     599          PVIEW(i) = 1.66
     600          PPSOL(i) = paprs(iof+i,1)
     601          zx_alpha1 = (paprs(iof+i,1)-pplay(iof+i,2))/(pplay(iof+i,1)-pplay(iof+i,2))
     602          zx_alpha2 = 1.0 - zx_alpha1
     603          PTL(i,1) = t(iof+i,1) * zx_alpha1 + t(iof+i,2) * zx_alpha2
     604          PTL(i,KLEV+1) = t(iof+i,KLEV)
     605          PDT0(i) = tsol(iof+i) - PTL(i,1)
     606       ENDDO
     607       DO k = 2, kflev
     608          DO i = 1, kdlon
     609             PTL(i,k) = (t(iof+i,k)+t(iof+i,k-1))*0.5
     610          ENDDO
     611       ENDDO
     612       DO k = 1, kflev
     613          DO i = 1, kdlon
     614             PDP(i,k) = paprs(iof+i,k)-paprs(iof+i,k+1)
     615             PTAVE(i,k) = t(iof+i,k)
     616             PWV(i,k) = MAX (q(iof+i,k), 1.0e-12)
     617             PQS(i,k) = PWV(i,k)
     618             !       Confert from  column density of ozone in a cell, in kDU, to a mass fraction
     619             POZON(i,k, :) = wo(iof+i, k, :) * RG * dobson_u * 1e3 &
     620                  / (paprs(iof+i, k) - paprs(iof+i, k+1))
     621             !       A activer pour CCMVAL on prend l'ozone impose (MPL 07042010)
     622             !       POZON(i,k,:) = wo(i,k,:) 
     623             !       print *,'RADLWSW: POZON',k, POZON(i,k,1)
     624             PCLDLD(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)
     625             PCLDLU(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)
     626             PCLDSW(i,k) = cldfra(iof+i,k)
     627             PTAU(i,1,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! 1e-12 serait instable
     628             PTAU(i,2,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! pour 32-bit machines
     629             POMEGA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAU(i,1,k))
     630             POMEGA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAU(i,2,k))
     631             PCG(i,1,k) = 0.865
     632             PCG(i,2,k) = 0.910
     633             !-
     634             ! Introduced for aerosol indirect forcings.
     635             ! The following values use the cloud optical thickness calculated from
     636             ! present-day aerosol concentrations whereas the quantities without the
     637             ! "A" at the end are for pre-industial (natural-only) aerosol concentrations
     638             !
     639             PTAUA(i,1,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! 1e-12 serait instable
     640             PTAUA(i,2,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! pour 32-bit machines
     641             POMEGAA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAUA(i,1,k))
     642             POMEGAA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAUA(i,2,k))
     643          ENDDO
     644       ENDDO
     645
     646       IF (type_trac == 'repr') THEN
     647#ifdef REPROBUS
     648          ndimozon = size(wo, 3)
     649          CALL RAD_INTERACTIF(POZON,iof)
     650#endif
     651       ENDIF
     652       !
     653       DO k = 1, kflev+1
     654          DO i = 1, kdlon
     655             PPMB(i,k) = paprs(iof+i,k)/100.0
     656          ENDDO
     657       ENDDO
     658       !
    659659!!!!! Modif MPL 6.01.09 avec RRTM, on passe de 5 a 6
    660     DO kk = 1, 6
    661       DO k = 1, kflev
    662         DO i = 1, kdlon
    663           PAER(i,k,kk) = 1.0E-15   !!!!! A REVOIR (MPL)
    664         ENDDO
    665       ENDDO
    666     ENDDO
    667     DO k = 1, kflev
    668       DO i = 1, kdlon
    669         tauaero(i,k,:,1)=tau_aero(iof+i,k,:,1)
    670         pizaero(i,k,:,1)=piz_aero(iof+i,k,:,1)
    671         cgaero(i,k,:,1) =cg_aero(iof+i,k,:,1)
    672         tauaero(i,k,:,2)=tau_aero(iof+i,k,:,2)
    673         pizaero(i,k,:,2)=piz_aero(iof+i,k,:,2)
    674         cgaero(i,k,:,2) =cg_aero(iof+i,k,:,2)
    675       ENDDO
    676     ENDDO
    677 !
    678 !===== iflag_rrtm ================================================
    679 !     
    680     IF (iflag_rrtm == 0) THEN       !!!! remettre 0 juste pour tester l'ancien rayt via rrtm
    681 !
    682 !--- Mise a zero des tableaux output du rayonnement LW-AR4 ----------             
    683       DO k = 1, kflev+1
    684       DO i = 1, kdlon
    685 !     print *,'RADLWSW: boucle mise a zero i k',i,k
    686       ZFLUP(i,k)=0.
    687       ZFLDN(i,k)=0.
    688       ZFLUP0(i,k)=0.
    689       ZFLDN0(i,k)=0.
    690       ZLWFT0_i(i,k)=0.
    691       ZFLUCUP_i(i,k)=0.
    692       ZFLUCDWN_i(i,k)=0.
    693       ENDDO
    694       ENDDO
    695       DO k = 1, kflev
    696          DO i = 1, kdlon
    697             zcool(i,k)=0.
    698             zcool_volc(i,k)=0. !NL
    699             zcool0(i,k)=0.
    700          ENDDO
    701       ENDDO
    702       DO i = 1, kdlon
    703       ztoplw(i)=0.
    704       zsollw(i)=0.
    705       ztoplw0(i)=0.
    706       zsollw0(i)=0.
    707       zsollwdown(i)=0.
    708       ENDDO
    709        ! Old radiation scheme, used for AR4 runs
    710        ! average day-night ozone for longwave
    711        CALL LW_LMDAR4(&
    712             PPMB, PDP,&
    713             PPSOL,PDT0,PEMIS,&
    714             PTL, PTAVE, PWV, POZON(:, :, 1), PAER,&
    715             PCLDLD,PCLDLU,&
    716             PVIEW,&
    717             zcool, zcool0,&
    718             ztoplw,zsollw,ztoplw0,zsollw0,&
    719             zsollwdown,&
    720             ZFLUP, ZFLDN, ZFLUP0,ZFLDN0)
    721 !----- Mise a zero des tableaux output du rayonnement SW-AR4
    722       DO k = 1, kflev+1
    723          DO i = 1, kdlon
    724             ZFSUP(i,k)=0.
    725             ZFSDN(i,k)=0.
    726             ZFSUP0(i,k)=0.
    727             ZFSDN0(i,k)=0.
    728             ZFSUPC0(i,k)=0.
    729             ZFSDNC0(i,k)=0.
    730             ZFLUPC0(i,k)=0.
    731             ZFLDNC0(i,k)=0.
    732             ZSWFT0_i(i,k)=0.
    733             ZFCUP_i(i,k)=0.
    734             ZFCDWN_i(i,k)=0.
    735             ZFCCUP_i(i,k)=0.
    736             ZFCCDWN_i(i,k)=0.
    737             ZFLCCUP_i(i,k)=0.
    738             ZFLCCDWN_i(i,k)=0.
    739             zswadaero(i,k)=0. !--NL
    740          ENDDO
    741       ENDDO
    742       DO k = 1, kflev
    743          DO i = 1, kdlon
    744             zheat(i,k)=0.
    745             zheat_volc(i,k)=0.
    746             zheat0(i,k)=0.
    747          ENDDO
    748       ENDDO
    749       DO i = 1, kdlon
    750       zalbpla(i)=0.
    751       ztopsw(i)=0.
    752       zsolsw(i)=0.
    753       ztopsw0(i)=0.
    754       zsolsw0(i)=0.
    755       ztopswadaero(i)=0.
    756       zsolswadaero(i)=0.
    757       ztopswaiaero(i)=0.
    758       zsolswaiaero(i)=0.
    759       ENDDO
    760 
    761       !--fraction of diffuse radiation in surface SW downward radiation
    762       !--not computed with old radiation scheme
    763       zsolswfdiff(:) = -999.999
    764 
    765 !     print *,'Avant SW_LMDAR4: PSCT zrmu0 zfract',PSCT, zrmu0, zfract
    766        ! daylight ozone, if we have it, for short wave
    767       CALL SW_AEROAR4(PSCT, zrmu0, zfract,&
     660       DO kk = 1, 6
     661          DO k = 1, kflev
     662             DO i = 1, kdlon
     663                PAER(i,k,kk) = 1.0E-15   !!!!! A REVOIR (MPL)
     664             ENDDO
     665          ENDDO
     666       ENDDO
     667       DO k = 1, kflev
     668          DO i = 1, kdlon
     669             tauaero(i,k,:,1)=tau_aero(iof+i,k,:,1)
     670             pizaero(i,k,:,1)=piz_aero(iof+i,k,:,1)
     671             cgaero(i,k,:,1) =cg_aero(iof+i,k,:,1)
     672             tauaero(i,k,:,2)=tau_aero(iof+i,k,:,2)
     673             pizaero(i,k,:,2)=piz_aero(iof+i,k,:,2)
     674             cgaero(i,k,:,2) =cg_aero(iof+i,k,:,2)
     675          ENDDO
     676       ENDDO
     677       !
     678       !===== iflag_rrtm ================================================
     679       !     
     680       IF (iflag_rrtm == 0) THEN       !!!! remettre 0 juste pour tester l'ancien rayt via rrtm
     681          !
     682          !--- Mise a zero des tableaux output du rayonnement LW-AR4 ----------             
     683          DO k = 1, kflev+1
     684             DO i = 1, kdlon
     685                !     print *,'RADLWSW: boucle mise a zero i k',i,k
     686                ZFLUP(i,k)=0.
     687                ZFLDN(i,k)=0.
     688                ZFLUP0(i,k)=0.
     689                ZFLDN0(i,k)=0.
     690                ZLWFT0_i(i,k)=0.
     691                ZFLUCUP_i(i,k)=0.
     692                ZFLUCDWN_i(i,k)=0.
     693             ENDDO
     694          ENDDO
     695          DO k = 1, kflev
     696             DO i = 1, kdlon
     697                zcool(i,k)=0.
     698                zcool_volc(i,k)=0. !NL
     699                zcool0(i,k)=0.
     700             ENDDO
     701          ENDDO
     702          DO i = 1, kdlon
     703             ztoplw(i)=0.
     704             zsollw(i)=0.
     705             ztoplw0(i)=0.
     706             zsollw0(i)=0.
     707             zsollwdown(i)=0.
     708          ENDDO
     709          ! Old radiation scheme, used for AR4 runs
     710          ! average day-night ozone for longwave
     711          CALL LW_LMDAR4(&
     712               PPMB, PDP,&
     713               PPSOL,PDT0,PEMIS,&
     714               PTL, PTAVE, PWV, POZON(:, :, 1), PAER,&
     715               PCLDLD,PCLDLU,&
     716               PVIEW,&
     717               zcool, zcool0,&
     718               ztoplw,zsollw,ztoplw0,zsollw0,&
     719               zsollwdown,&
     720               ZFLUP, ZFLDN, ZFLUP0,ZFLDN0)
     721          !----- Mise a zero des tableaux output du rayonnement SW-AR4
     722          DO k = 1, kflev+1
     723             DO i = 1, kdlon
     724                ZFSUP(i,k)=0.
     725                ZFSDN(i,k)=0.
     726                ZFSUP0(i,k)=0.
     727                ZFSDN0(i,k)=0.
     728                ZFSUPC0(i,k)=0.
     729                ZFSDNC0(i,k)=0.
     730                ZFLUPC0(i,k)=0.
     731                ZFLDNC0(i,k)=0.
     732                ZSWFT0_i(i,k)=0.
     733                ZFCUP_i(i,k)=0.
     734                ZFCDWN_i(i,k)=0.
     735                ZFCCUP_i(i,k)=0.
     736                ZFCCDWN_i(i,k)=0.
     737                ZFLCCUP_i(i,k)=0.
     738                ZFLCCDWN_i(i,k)=0.
     739                zswadaero(i,k)=0. !--NL
     740             ENDDO
     741          ENDDO
     742          DO k = 1, kflev
     743             DO i = 1, kdlon
     744                zheat(i,k)=0.
     745                zheat_volc(i,k)=0.
     746                zheat0(i,k)=0.
     747             ENDDO
     748          ENDDO
     749          DO i = 1, kdlon
     750             zalbpla(i)=0.
     751             ztopsw(i)=0.
     752             zsolsw(i)=0.
     753             ztopsw0(i)=0.
     754             zsolsw0(i)=0.
     755             ztopswadaero(i)=0.
     756             zsolswadaero(i)=0.
     757             ztopswaiaero(i)=0.
     758             zsolswaiaero(i)=0.
     759          ENDDO
     760
     761          !--fraction of diffuse radiation in surface SW downward radiation
     762          !--not computed with old radiation scheme
     763          zsolswfdiff(:) = -999.999
     764
     765          !     print *,'Avant SW_LMDAR4: PSCT zrmu0 zfract',PSCT, zrmu0, zfract
     766          ! daylight ozone, if we have it, for short wave
     767          CALL SW_AEROAR4(PSCT, zrmu0, zfract,&
    768768               PPMB, PDP,&
    769769               PPSOL, PALBD, PALBP,&
     
    783783               ok_ade, ok_aie, flag_aerosol,flag_aerosol_strat)
    784784
    785        ZSWFT0_i(:,:) = ZFSDN0(:,:)-ZFSUP0(:,:)
    786        ZLWFT0_i(:,:) =-ZFLDN0(:,:)-ZFLUP0(:,:)
    787 
    788        DO i=1,kdlon
    789        DO k=1,kflev+1
    790          lwdn0 ( iof+i,k)   = ZFLDN0 ( i,k)
    791          lwdn  ( iof+i,k)   = ZFLDN  ( i,k)
    792          lwup0 ( iof+i,k)   = ZFLUP0 ( i,k)
    793          lwup  ( iof+i,k)   = ZFLUP  ( i,k)
    794          swdn0 ( iof+i,k)   = ZFSDN0 ( i,k)
    795          swdn  ( iof+i,k)   = ZFSDN  ( i,k)
    796          swup0 ( iof+i,k)   = ZFSUP0 ( i,k)
    797          swup  ( iof+i,k)   = ZFSUP  ( i,k)
    798        ENDDO 
    799        ENDDO 
    800 !
    801     ELSE IF (iflag_rrtm == 1) then 
     785          ZSWFT0_i(:,:) = ZFSDN0(:,:)-ZFSUP0(:,:)
     786          ZLWFT0_i(:,:) =-ZFLDN0(:,:)-ZFLUP0(:,:)
     787
     788          DO i=1,kdlon
     789             DO k=1,kflev+1
     790                lwdn0 ( iof+i,k)   = ZFLDN0 ( i,k)
     791                lwdn  ( iof+i,k)   = ZFLDN  ( i,k)
     792                lwup0 ( iof+i,k)   = ZFLUP0 ( i,k)
     793                lwup  ( iof+i,k)   = ZFLUP  ( i,k)
     794                swdn0 ( iof+i,k)   = ZFSDN0 ( i,k)
     795                swdn  ( iof+i,k)   = ZFSDN  ( i,k)
     796                swup0 ( iof+i,k)   = ZFSUP0 ( i,k)
     797                swup  ( iof+i,k)   = ZFSUP  ( i,k)
     798             ENDDO
     799          ENDDO
     800          !
     801       ELSE IF (iflag_rrtm == 1) then 
    802802#ifdef CPP_RRTM
    803 !      if (prt_level.gt.10)write(lunout,*)'CPP_RRTM=.T.'
    804 !===== iflag_rrtm=1, on passe dans SW via RECMWFL ===============
    805 
    806       DO k = 1, kflev+1
    807       DO i = 1, kdlon
    808         ZEMTD_i(i,k)=0.
    809         ZEMTU_i(i,k)=0.
    810         ZTRSO_i(i,k)=0.
    811         ZTH_i(i,k)=0.
    812         ZLWFT_i(i,k)=0.
    813         ZSWFT_i(i,k)=0.
    814         ZFLUX_i(i,1,k)=0.
    815         ZFLUX_i(i,2,k)=0.
    816         ZFLUC_i(i,1,k)=0.
    817         ZFLUC_i(i,2,k)=0.
    818         ZFSDWN_i(i,k)=0.
    819         ZFCDWN_i(i,k)=0.
    820         ZFCCDWN_i(i,k)=0.
    821         ZFSUP_i(i,k)=0.
    822         ZFCUP_i(i,k)=0.
    823         ZFCCUP_i(i,k)=0.
    824         ZFLCCDWN_i(i,k)=0.
    825         ZFLCCUP_i(i,k)=0.
    826       ENDDO
    827       ENDDO
    828 !
    829 !--OB
    830 !--aerosol TOT  - anthropogenic+natural - index 2
    831 !--aerosol NAT  - natural only          - index 1
    832 !
    833       DO i = 1, kdlon
    834       DO k = 1, kflev
    835       DO kk=1, NSW
    836 !
    837       PTAU_TOT(i,kflev+1-k,kk)=tau_aero_sw_rrtm(i,k,2,kk)
    838       PPIZA_TOT(i,kflev+1-k,kk)=piz_aero_sw_rrtm(i,k,2,kk)
    839       PCGA_TOT(i,kflev+1-k,kk)=cg_aero_sw_rrtm(i,k,2,kk)
    840 !
    841       PTAU_NAT(i,kflev+1-k,kk)=tau_aero_sw_rrtm(i,k,1,kk)
    842       PPIZA_NAT(i,kflev+1-k,kk)=piz_aero_sw_rrtm(i,k,1,kk)
    843       PCGA_NAT(i,kflev+1-k,kk)=cg_aero_sw_rrtm(i,k,1,kk)
    844 !
    845       ENDDO
    846       ENDDO
    847       ENDDO
    848 !-end OB
    849 !
    850 !--C. Kleinschmitt
    851 !--aerosol TOT  - anthropogenic+natural - index 2
    852 !--aerosol NAT  - natural only          - index 1
    853 !
    854       DO i = 1, kdlon
    855       DO k = 1, kflev
    856       DO kk=1, NLW
    857 !
    858       PTAU_LW_TOT(i,kflev+1-k,kk)=tau_aero_lw_rrtm(i,k,2,kk)
    859       PTAU_LW_NAT(i,kflev+1-k,kk)=tau_aero_lw_rrtm(i,k,1,kk)
    860 !
    861       ENDDO
    862       ENDDO
    863       ENDDO
    864 !-end C. Kleinschmitt
    865 !     
    866       DO i = 1, kdlon
    867       ZCTRSO(i,1)=0.
    868       ZCTRSO(i,2)=0.
    869       ZCEMTR(i,1)=0.
    870       ZCEMTR(i,2)=0.
    871       ZTRSOD(i)=0.
    872       ZLWFC(i,1)=0.
    873       ZLWFC(i,2)=0.
    874       ZSWFC(i,1)=0.
    875       ZSWFC(i,2)=0.
    876       PFSDNN(i)=0.
    877       PFSDNV(i)=0.
    878       DO kk = 1, NSW
    879         PSFSWDIR(i,kk)=0.
    880         PSFSWDIF(i,kk)=0.
    881       ENDDO
    882       ENDDO
    883 !----- Fin des mises a zero des tableaux output de RECMWF -------------------             
    884 !        GEMU(1:klon)=sin(rlatd(1:klon))
    885 ! On met les donnees dans l'ordre des niveaux arpege
    886          paprs_i(:,1)=paprs(:,klev+1)
    887          DO k=1,klev
    888             paprs_i(1:klon,k+1) =paprs(1:klon,klev+1-k)
    889             pplay_i(1:klon,k)   =pplay(1:klon,klev+1-k)
    890             cldfra_i(1:klon,k)  =cldfra(1:klon,klev+1-k)
    891             PDP_i(1:klon,k)     =PDP(1:klon,klev+1-k)
    892             t_i(1:klon,k)       =t(1:klon,klev+1-k)
    893             q_i(1:klon,k)       =q(1:klon,klev+1-k)
    894             qsat_i(1:klon,k)    =qsat(1:klon,klev+1-k)
    895             flwc_i(1:klon,k)    =flwc(1:klon,klev+1-k)
    896             fiwc_i(1:klon,k)    =fiwc(1:klon,klev+1-k)
    897             ref_liq_i(1:klon,k) =ref_liq(1:klon,klev+1-k)
    898             ref_ice_i(1:klon,k) =ref_ice(1:klon,klev+1-k)
    899 !-OB
    900             ref_liq_pi_i(1:klon,k) =ref_liq_pi(1:klon,klev+1-k)
    901             ref_ice_pi_i(1:klon,k) =ref_ice_pi(1:klon,klev+1-k)
    902          ENDDO
    903          DO k=1,kflev
    904            POZON_i(1:klon,k,:)=POZON(1:klon,kflev+1-k,:)
     803          !      if (prt_level.gt.10)write(lunout,*)'CPP_RRTM=.T.'
     804          !===== iflag_rrtm=1, on passe dans SW via RECMWFL ===============
     805
     806          DO k = 1, kflev+1
     807             DO i = 1, kdlon
     808                ZEMTD_i(i,k)=0.
     809                ZEMTU_i(i,k)=0.
     810                ZTRSO_i(i,k)=0.
     811                ZTH_i(i,k)=0.
     812                ZLWFT_i(i,k)=0.
     813                ZSWFT_i(i,k)=0.
     814                ZFLUX_i(i,1,k)=0.
     815                ZFLUX_i(i,2,k)=0.
     816                ZFLUC_i(i,1,k)=0.
     817                ZFLUC_i(i,2,k)=0.
     818                ZFSDWN_i(i,k)=0.
     819                ZFCDWN_i(i,k)=0.
     820                ZFCCDWN_i(i,k)=0.
     821                ZFSUP_i(i,k)=0.
     822                ZFCUP_i(i,k)=0.
     823                ZFCCUP_i(i,k)=0.
     824                ZFLCCDWN_i(i,k)=0.
     825                ZFLCCUP_i(i,k)=0.
     826             ENDDO
     827          ENDDO
     828          !
     829          !--OB
     830          !--aerosol TOT  - anthropogenic+natural - index 2
     831          !--aerosol NAT  - natural only          - index 1
     832          !
     833          DO i = 1, kdlon
     834             DO k = 1, kflev
     835                DO kk=1, NSW
     836                   !
     837                   PTAU_TOT(i,kflev+1-k,kk)=tau_aero_sw_rrtm(i,k,2,kk)
     838                   PPIZA_TOT(i,kflev+1-k,kk)=piz_aero_sw_rrtm(i,k,2,kk)
     839                   PCGA_TOT(i,kflev+1-k,kk)=cg_aero_sw_rrtm(i,k,2,kk)
     840                   !
     841                   PTAU_NAT(i,kflev+1-k,kk)=tau_aero_sw_rrtm(i,k,1,kk)
     842                   PPIZA_NAT(i,kflev+1-k,kk)=piz_aero_sw_rrtm(i,k,1,kk)
     843                   PCGA_NAT(i,kflev+1-k,kk)=cg_aero_sw_rrtm(i,k,1,kk)
     844                   !
     845                ENDDO
     846             ENDDO
     847          ENDDO
     848          !-end OB
     849          !
     850          !--C. Kleinschmitt
     851          !--aerosol TOT  - anthropogenic+natural - index 2
     852          !--aerosol NAT  - natural only          - index 1
     853          !
     854          DO i = 1, kdlon
     855             DO k = 1, kflev
     856                DO kk=1, NLW
     857                   !
     858                   PTAU_LW_TOT(i,kflev+1-k,kk)=tau_aero_lw_rrtm(i,k,2,kk)
     859                   PTAU_LW_NAT(i,kflev+1-k,kk)=tau_aero_lw_rrtm(i,k,1,kk)
     860                   !
     861                ENDDO
     862             ENDDO
     863          ENDDO
     864          !-end C. Kleinschmitt
     865          !     
     866          DO i = 1, kdlon
     867             ZCTRSO(i,1)=0.
     868             ZCTRSO(i,2)=0.
     869             ZCEMTR(i,1)=0.
     870             ZCEMTR(i,2)=0.
     871             ZTRSOD(i)=0.
     872             ZLWFC(i,1)=0.
     873             ZLWFC(i,2)=0.
     874             ZSWFC(i,1)=0.
     875             ZSWFC(i,2)=0.
     876             PFSDNN(i)=0.
     877             PFSDNV(i)=0.
     878             DO kk = 1, NSW
     879                PSFSWDIR(i,kk)=0.
     880                PSFSWDIF(i,kk)=0.
     881             ENDDO
     882          ENDDO
     883          !----- Fin des mises a zero des tableaux output de RECMWF -------------------             
     884          !        GEMU(1:klon)=sin(rlatd(1:klon))
     885          ! On met les donnees dans l'ordre des niveaux arpege
     886          paprs_i(:,1)=paprs(:,klev+1)
     887          DO k=1,klev
     888             paprs_i(1:klon,k+1) =paprs(1:klon,klev+1-k)
     889             pplay_i(1:klon,k)   =pplay(1:klon,klev+1-k)
     890             cldfra_i(1:klon,k)  =cldfra(1:klon,klev+1-k)
     891             PDP_i(1:klon,k)     =PDP(1:klon,klev+1-k)
     892             t_i(1:klon,k)       =t(1:klon,klev+1-k)
     893             q_i(1:klon,k)       =q(1:klon,klev+1-k)
     894             qsat_i(1:klon,k)    =qsat(1:klon,klev+1-k)
     895             flwc_i(1:klon,k)    =flwc(1:klon,klev+1-k)
     896             fiwc_i(1:klon,k)    =fiwc(1:klon,klev+1-k)
     897             ref_liq_i(1:klon,k) =ref_liq(1:klon,klev+1-k)
     898             ref_ice_i(1:klon,k) =ref_ice(1:klon,klev+1-k)
     899             !-OB
     900             ref_liq_pi_i(1:klon,k) =ref_liq_pi(1:klon,klev+1-k)
     901             ref_ice_pi_i(1:klon,k) =ref_ice_pi(1:klon,klev+1-k)
     902          ENDDO
     903          DO k=1,kflev
     904             POZON_i(1:klon,k,:)=POZON(1:klon,kflev+1-k,:)
    905905!!!            POZON_i(1:klon,k)=POZON(1:klon,k)            !!! on laisse 1=sol et klev=top
    906 !          print *,'Juste avant RECMWFL: k tsol temp',k,tsol,t(1,k)
     906             !          print *,'Juste avant RECMWFL: k tsol temp',k,tsol,t(1,k)
    907907!!!!!!! Modif MPL 6.01.09 avec RRTM, on passe de 5 a 6     
    908             DO i=1,6
    909             PAER_i(1:klon,k,i)=PAER(1:klon,kflev+1-k,i)
    910             ENDDO
    911          ENDDO
    912 
    913 !       print *,'RADLWSW: avant RECMWFL, RI0,rmu0=',solaire,rmu0
    914 
    915 !  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    916 ! La version ARPEGE1D utilise differentes valeurs de la constante
    917 ! solaire suivant le rayonnement utilise.
    918 ! A controler ...
    919 ! SOLAR FLUX AT THE TOP (/YOMPHY3/)
    920 ! introduce season correction
    921 !--------------------------------------
    922 ! RII0 = RIP0
    923 ! IF(LRAYFM)
    924 ! RII0 = RIP0M   ! =rip0m if Morcrette non-each time step call.
    925 ! IF(LRAYFM15)
    926 ! RII0 = RIP0M15 ! =rip0m if Morcrette non-each time step call.
    927          RII0=solaire/zdist/zdist
    928 !  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    929 ! Ancien appel a RECMWF (celui du cy25)
    930 !        CALL RECMWF (ist , iend, klon , ktdia , klev   , kmode ,
    931 !    s   PALBD    , PALBP   , paprs_i , pplay_i , RCO2   , cldfra_i,
    932 !    s   POZON_i  , PAER_i  , PDP_i   , PEMIS   , GEMU   , rmu0,
    933 !    s    q_i     , qsat_i  , fiwc_i  , flwc_i  , zmasq  , t_i  ,tsol,
    934 !    s   ZEMTD_i  , ZEMTU_i , ZTRSO_i ,
    935 !    s   ZTH_i    , ZCTRSO  , ZCEMTR  , ZTRSOD  ,
    936 !    s   ZLWFC    , ZLWFT_i , ZSWFC   , ZSWFT_i ,
    937 !    s   ZFLUX_i  , ZFLUC_i , ZFSDWN_i, ZFSUP_i , ZFCDWN_i,ZFCUP_i)
    938 !    s   'RECMWF ')
    939 !
    940       IF (lldebug) THEN
    941         CALL writefield_phy('paprs_i',paprs_i,klev+1)
    942         CALL writefield_phy('pplay_i',pplay_i,klev)
    943         CALL writefield_phy('cldfra_i',cldfra_i,klev)
    944         CALL writefield_phy('pozon_i',POZON_i,klev)
    945         CALL writefield_phy('paer_i',PAER_i,klev)
    946         CALL writefield_phy('pdp_i',PDP_i,klev)
    947         CALL writefield_phy('q_i',q_i,klev)
    948         CALL writefield_phy('qsat_i',qsat_i,klev)
    949         CALL writefield_phy('fiwc_i',fiwc_i,klev)
    950         CALL writefield_phy('flwc_i',flwc_i,klev)
    951         CALL writefield_phy('t_i',t_i,klev)
    952         CALL writefield_phy('palbd_new',PALBD_NEW,NSW)
    953         CALL writefield_phy('palbp_new',PALBP_NEW,NSW)
    954       ENDIF
    955 
    956 ! Nouvel appel a RECMWF (celui du cy32t0)
    957          CALL RECMWF_AERO (ist , iend, klon , ktdia  , klev   , kmode ,&
    958          PALBD_NEW,PALBP_NEW, paprs_i , pplay_i , RCO2   , cldfra_i,&
    959          POZON_i  , PAER_i  , PDP_i   , PEMIS   , rmu0   ,&
    960          q_i     , qsat_i  , fiwc_i  , flwc_i  , zmasq  , t_i  ,tsol,&
    961          ref_liq_i, ref_ice_i, &
    962          ref_liq_pi_i, ref_ice_pi_i, &   ! rajoute par OB pour diagnostiquer effet indirect
    963          ZEMTD_i  , ZEMTU_i , ZTRSO_i ,&
    964          ZTH_i    , ZCTRSO  , ZCEMTR  , ZTRSOD  ,&
    965          ZLWFC    , ZLWFT_i , ZSWFC   , ZSWFT_i ,&
    966          PSFSWDIR , PSFSWDIF, PFSDNN  , PFSDNV  ,&
    967          PPIZA_TOT, PCGA_TOT,PTAU_TOT,&
    968          PPIZA_NAT, PCGA_NAT,PTAU_NAT,           &  ! rajoute par OB pour diagnostiquer effet direct
    969          PTAU_LW_TOT, PTAU_LW_NAT,               &  ! rajoute par C. Kleinschmitt
    970          ZFLUX_i  , ZFLUC_i ,&
    971          ZFSDWN_i , ZFSUP_i , ZFCDWN_i, ZFCUP_i, ZFCCDWN_i, ZFCCUP_i, ZFLCCDWN_i, ZFLCCUP_i, &
    972          ZTOPSWADAERO,ZSOLSWADAERO,&  ! rajoute par OB pour diagnostics
    973          ZTOPSWAD0AERO,ZSOLSWAD0AERO,&
    974          ZTOPSWAIAERO,ZSOLSWAIAERO, &
    975          ZTOPSWCF_AERO,ZSOLSWCF_AERO, &
    976          ZSWADAERO, & !--NL
    977          ZTOPLWADAERO,ZSOLLWADAERO,&  ! rajoute par C. Kleinscmitt pour LW diagnostics
    978          ZTOPLWAD0AERO,ZSOLLWAD0AERO,&
    979          ZTOPLWAIAERO,ZSOLLWAIAERO, &
    980          ZLWADAERO, & !--NL
    981          volmip_solsw, flag_volc_surfstrat, & !--VOLMIP
    982          ok_ade, ok_aie, ok_volcan, flag_aerosol,flag_aerosol_strat, flag_aer_feedback) ! flags aerosols
    983 
    984 !--OB diagnostics
    985 ! & PTOPSWAIAERO,PSOLSWAIAERO,&
    986 ! & PTOPSWCFAERO,PSOLSWCFAERO,&
    987 ! & PSWADAERO,& !--NL
    988 !!--LW diagnostics CK
    989 ! & PTOPLWADAERO,PSOLLWADAERO,&
    990 ! & PTOPLWAD0AERO,PSOLLWAD0AERO,&
    991 ! & PTOPLWAIAERO,PSOLLWAIAERO,&
    992 ! & PLWADAERO,& !--NL
    993 !!..end
    994 ! & ok_ade, ok_aie, ok_volcan, flag_aerosol,flag_aerosol_strat,&
    995 ! & flag_aer_feedback)
    996 
    997            
    998 !        print *,'RADLWSW: apres RECMWF'
    999       IF (lldebug) THEN
    1000         CALL writefield_phy('zemtd_i',ZEMTD_i,klev+1)
    1001         CALL writefield_phy('zemtu_i',ZEMTU_i,klev+1)
    1002         CALL writefield_phy('ztrso_i',ZTRSO_i,klev+1)
    1003         CALL writefield_phy('zth_i',ZTH_i,klev+1)
    1004         CALL writefield_phy('zctrso',ZCTRSO,2)
    1005         CALL writefield_phy('zcemtr',ZCEMTR,2)
    1006         CALL writefield_phy('ztrsod',ZTRSOD,1)
    1007         CALL writefield_phy('zlwfc',ZLWFC,2)
    1008         CALL writefield_phy('zlwft_i',ZLWFT_i,klev+1)
    1009         CALL writefield_phy('zswfc',ZSWFC,2)
    1010         CALL writefield_phy('zswft_i',ZSWFT_i,klev+1)
    1011         CALL writefield_phy('psfswdir',PSFSWDIR,6)
    1012         CALL writefield_phy('psfswdif',PSFSWDIF,6)
    1013         CALL writefield_phy('pfsdnn',PFSDNN,1)
    1014         CALL writefield_phy('pfsdnv',PFSDNV,1)
    1015         CALL writefield_phy('ppiza_dst',PPIZA_TOT,klev)
    1016         CALL writefield_phy('pcga_dst',PCGA_TOT,klev)
    1017         CALL writefield_phy('ptaurel_dst',PTAU_TOT,klev)
    1018         CALL writefield_phy('zflux_i',ZFLUX_i,klev+1)
    1019         CALL writefield_phy('zfluc_i',ZFLUC_i,klev+1)
    1020         CALL writefield_phy('zfsdwn_i',ZFSDWN_i,klev+1)
    1021         CALL writefield_phy('zfsup_i',ZFSUP_i,klev+1)
    1022         CALL writefield_phy('zfcdwn_i',ZFCDWN_i,klev+1)
    1023         CALL writefield_phy('zfcup_i',ZFCUP_i,klev+1)
    1024       ENDIF
    1025 
    1026 ! ---------
    1027 ! ---------
    1028 ! On retablit l'ordre des niveaux lmd pour les tableaux de sortie
    1029 ! D autre part, on multiplie les resultats SW par fract pour etre coherent
    1030 ! avec l ancien rayonnement AR4. Si nuit, fract=0 donc pas de
    1031 ! rayonnement SW. (MPL 260609)
    1032       DO k=0,klev
    1033          DO i=1,klon
    1034          ZEMTD(i,k+1)  = ZEMTD_i(i,k+1)
    1035          ZEMTU(i,k+1)  = ZEMTU_i(i,k+1)
    1036          ZTRSO(i,k+1)  = ZTRSO_i(i,k+1)
    1037          ZTH(i,k+1)    = ZTH_i(i,k+1)
    1038 !        ZLWFT(i,k+1)  = ZLWFT_i(i,klev+1-k)
    1039 !        ZSWFT(i,k+1)  = ZSWFT_i(i,klev+1-k)
    1040          ZFLUP(i,k+1)  = ZFLUX_i(i,1,k+1)
    1041          ZFLDN(i,k+1)  = ZFLUX_i(i,2,k+1)
    1042          ZFLUP0(i,k+1) = ZFLUC_i(i,1,k+1)
    1043          ZFLDN0(i,k+1) = ZFLUC_i(i,2,k+1)
    1044          ZFSDN(i,k+1)  = ZFSDWN_i(i,k+1)*fract(i)
    1045          ZFSDN0(i,k+1) = ZFCDWN_i(i,k+1)*fract(i)
    1046          ZFSDNC0(i,k+1)= ZFCCDWN_i(i,k+1)*fract(i)
    1047          ZFSUP (i,k+1) = ZFSUP_i(i,k+1)*fract(i)
    1048          ZFSUP0(i,k+1) = ZFCUP_i(i,k+1)*fract(i)
    1049          ZFSUPC0(i,k+1)= ZFCCUP_i(i,k+1)*fract(i)
    1050          ZFLDNC0(i,k+1)= ZFLCCDWN_i(i,k+1)
    1051          ZFLUPC0(i,k+1)= ZFLCCUP_i(i,k+1)
    1052          IF (ok_volcan) THEN
    1053             ZSWADAERO(i,k+1)=ZSWADAERO(i,k+1)*fract(i) !--NL
    1054          ENDIF
    1055          
    1056 !   Nouveau calcul car visiblement ZSWFT et ZSWFC sont nuls dans RRTM cy32
    1057 !   en sortie de radlsw.F90 - MPL 7.01.09
    1058          ZSWFT(i,k+1)  = (ZFSDWN_i(i,k+1)-ZFSUP_i(i,k+1))*fract(i)
    1059          ZSWFT0_i(i,k+1) = (ZFCDWN_i(i,k+1)-ZFCUP_i(i,k+1))*fract(i)
    1060 !        WRITE(*,'("FSDN FSUP FCDN FCUP: ",4E12.5)') ZFSDWN_i(i,k+1),&
    1061 !        ZFSUP_i(i,k+1),ZFCDWN_i(i,k+1),ZFCUP_i(i,k+1)
    1062          ZLWFT(i,k+1) =-ZFLUX_i(i,2,k+1)-ZFLUX_i(i,1,k+1)
    1063          ZLWFT0_i(i,k+1)=-ZFLUC_i(i,2,k+1)-ZFLUC_i(i,1,k+1)
    1064 !        print *,'FLUX2 FLUX1 FLUC2 FLUC1',ZFLUX_i(i,2,k+1),&
    1065 !    & ZFLUX_i(i,1,k+1),ZFLUC_i(i,2,k+1),ZFLUC_i(i,1,k+1)
    1066          ENDDO
    1067       ENDDO
    1068 
    1069 !--ajout OB
    1070       ZTOPSWADAERO(:) =ZTOPSWADAERO(:) *fract(:)
    1071       ZSOLSWADAERO(:) =ZSOLSWADAERO(:) *fract(:)
    1072       ZTOPSWAD0AERO(:)=ZTOPSWAD0AERO(:)*fract(:)
    1073       ZSOLSWAD0AERO(:)=ZSOLSWAD0AERO(:)*fract(:)
    1074       ZTOPSWAIAERO(:) =ZTOPSWAIAERO(:) *fract(:)
    1075       ZSOLSWAIAERO(:) =ZSOLSWAIAERO(:) *fract(:)
    1076       ZTOPSWCF_AERO(:,1)=ZTOPSWCF_AERO(:,1)*fract(:)
    1077       ZTOPSWCF_AERO(:,2)=ZTOPSWCF_AERO(:,2)*fract(:)
    1078       ZTOPSWCF_AERO(:,3)=ZTOPSWCF_AERO(:,3)*fract(:)
    1079       ZSOLSWCF_AERO(:,1)=ZSOLSWCF_AERO(:,1)*fract(:)
    1080       ZSOLSWCF_AERO(:,2)=ZSOLSWCF_AERO(:,2)*fract(:)
    1081       ZSOLSWCF_AERO(:,3)=ZSOLSWCF_AERO(:,3)*fract(:)
    1082 
    1083 ! ---------
    1084 ! ---------
    1085 ! On renseigne les champs LMDz, pour avoir la meme chose qu'en sortie de
    1086 ! LW_LMDAR4 et SW_LMDAR4
    1087 
    1088       !--fraction of diffuse radiation in surface SW downward radiation
    1089       DO i = 1, kdlon
    1090        IF (fract(i).GT.0.0) THEN
    1091          zdir=SUM(PSFSWDIR(i,:))
    1092          zdif=SUM(PSFSWDIF(i,:))
    1093          zsolswfdiff(i) = zdif/(zdir+zdif)
    1094        ELSE  !--night
    1095          zsolswfdiff(i) = 1.0
    1096        ENDIF
    1097       ENDDO
    1098 !
    1099       DO i = 1, kdlon
    1100          zsolsw(i)    = ZSWFT(i,1)
    1101          zsolsw0(i)   = ZSWFT0_i(i,1)
    1102 !        zsolsw0(i)   = ZFSDN0(i,1)     -ZFSUP0(i,1)
    1103          ztopsw(i)    = ZSWFT(i,klev+1)
    1104          ztopsw0(i)   = ZSWFT0_i(i,klev+1)
    1105 !        ztopsw0(i)   = ZFSDN0(i,klev+1)-ZFSUP0(i,klev+1)
    1106 !         
    1107 !        zsollw(i)    = ZFLDN(i,1)      -ZFLUP(i,1)
    1108 !        zsollw0(i)   = ZFLDN0(i,1)     -ZFLUP0(i,1)
    1109 !        ztoplw(i)    = ZFLDN(i,klev+1) -ZFLUP(i,klev+1)
    1110 !        ztoplw0(i)   = ZFLDN0(i,klev+1)-ZFLUP0(i,klev+1)
    1111          zsollw(i)    = ZLWFT(i,1)
    1112          zsollw0(i)   = ZLWFT0_i(i,1)
    1113          ztoplw(i)    = ZLWFT(i,klev+1)*(-1)
    1114          ztoplw0(i)   = ZLWFT0_i(i,klev+1)*(-1)
    1115 !         
    1116          IF (fract(i) == 0.) THEN
     908             DO i=1,6
     909                PAER_i(1:klon,k,i)=PAER(1:klon,kflev+1-k,i)
     910             ENDDO
     911          ENDDO
     912
     913          !       print *,'RADLWSW: avant RECMWFL, RI0,rmu0=',solaire,rmu0
     914
     915          !  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     916          ! La version ARPEGE1D utilise differentes valeurs de la constante
     917          ! solaire suivant le rayonnement utilise.
     918          ! A controler ...
     919          ! SOLAR FLUX AT THE TOP (/YOMPHY3/)
     920          ! introduce season correction
     921          !--------------------------------------
     922          ! RII0 = RIP0
     923          ! IF(LRAYFM)
     924          ! RII0 = RIP0M   ! =rip0m if Morcrette non-each time step call.
     925          ! IF(LRAYFM15)
     926          ! RII0 = RIP0M15 ! =rip0m if Morcrette non-each time step call.
     927          RII0=solaire/zdist/zdist
     928          !  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     929          ! Ancien appel a RECMWF (celui du cy25)
     930          !        CALL RECMWF (ist , iend, klon , ktdia , klev   , kmode ,
     931          !    s   PALBD    , PALBP   , paprs_i , pplay_i , RCO2   , cldfra_i,
     932          !    s   POZON_i  , PAER_i  , PDP_i   , PEMIS   , GEMU   , rmu0,
     933          !    s    q_i     , qsat_i  , fiwc_i  , flwc_i  , zmasq  , t_i  ,tsol,
     934          !    s   ZEMTD_i  , ZEMTU_i , ZTRSO_i ,
     935          !    s   ZTH_i    , ZCTRSO  , ZCEMTR  , ZTRSOD  ,
     936          !    s   ZLWFC    , ZLWFT_i , ZSWFC   , ZSWFT_i ,
     937          !    s   ZFLUX_i  , ZFLUC_i , ZFSDWN_i, ZFSUP_i , ZFCDWN_i,ZFCUP_i)
     938          !    s   'RECMWF ')
     939          !
     940          IF (lldebug) THEN
     941             CALL writefield_phy('paprs_i',paprs_i,klev+1)
     942             CALL writefield_phy('pplay_i',pplay_i,klev)
     943             CALL writefield_phy('cldfra_i',cldfra_i,klev)
     944             CALL writefield_phy('pozon_i',POZON_i,klev)
     945             CALL writefield_phy('paer_i',PAER_i,klev)
     946             CALL writefield_phy('pdp_i',PDP_i,klev)
     947             CALL writefield_phy('q_i',q_i,klev)
     948             CALL writefield_phy('qsat_i',qsat_i,klev)
     949             CALL writefield_phy('fiwc_i',fiwc_i,klev)
     950             CALL writefield_phy('flwc_i',flwc_i,klev)
     951             CALL writefield_phy('t_i',t_i,klev)
     952             CALL writefield_phy('palbd_new',PALBD_NEW,NSW)
     953             CALL writefield_phy('palbp_new',PALBP_NEW,NSW)
     954          ENDIF
     955
     956          ! Nouvel appel a RECMWF (celui du cy32t0)
     957          CALL RECMWF_AERO (ist , iend, klon , ktdia  , klev   , kmode ,&
     958               PALBD_NEW,PALBP_NEW, paprs_i , pplay_i , RCO2   , cldfra_i,&
     959               POZON_i  , PAER_i  , PDP_i   , PEMIS   , rmu0   ,&
     960               q_i     , qsat_i  , fiwc_i  , flwc_i  , zmasq  , t_i  ,tsol,&
     961               ref_liq_i, ref_ice_i, &
     962               ref_liq_pi_i, ref_ice_pi_i, &   ! rajoute par OB pour diagnostiquer effet indirect
     963               ZEMTD_i  , ZEMTU_i , ZTRSO_i ,&
     964               ZTH_i    , ZCTRSO  , ZCEMTR  , ZTRSOD  ,&
     965               ZLWFC    , ZLWFT_i , ZSWFC   , ZSWFT_i ,&
     966               PSFSWDIR , PSFSWDIF, PFSDNN  , PFSDNV  ,&
     967               PPIZA_TOT, PCGA_TOT,PTAU_TOT,&
     968               PPIZA_NAT, PCGA_NAT,PTAU_NAT,           &  ! rajoute par OB pour diagnostiquer effet direct
     969               PTAU_LW_TOT, PTAU_LW_NAT,               &  ! rajoute par C. Kleinschmitt
     970               ZFLUX_i  , ZFLUC_i ,&
     971               ZFSDWN_i , ZFSUP_i , ZFCDWN_i, ZFCUP_i, ZFCCDWN_i, ZFCCUP_i, ZFLCCDWN_i, ZFLCCUP_i, &
     972               ZTOPSWADAERO,ZSOLSWADAERO,&  ! rajoute par OB pour diagnostics
     973               ZTOPSWAD0AERO,ZSOLSWAD0AERO,&
     974               ZTOPSWAIAERO,ZSOLSWAIAERO, &
     975               ZTOPSWCF_AERO,ZSOLSWCF_AERO, &
     976               ZSWADAERO, & !--NL
     977               ZTOPLWADAERO,ZSOLLWADAERO,&  ! rajoute par C. Kleinscmitt pour LW diagnostics
     978               ZTOPLWAD0AERO,ZSOLLWAD0AERO,&
     979               ZTOPLWAIAERO,ZSOLLWAIAERO, &
     980               ZLWADAERO, & !--NL
     981               volmip_solsw, flag_volc_surfstrat, & !--VOLMIP
     982               ok_ade, ok_aie, ok_volcan, flag_aerosol,flag_aerosol_strat, flag_aer_feedback) ! flags aerosols
     983
     984          !--OB diagnostics
     985          ! & PTOPSWAIAERO,PSOLSWAIAERO,&
     986          ! & PTOPSWCFAERO,PSOLSWCFAERO,&
     987          ! & PSWADAERO,& !--NL
     988          !!--LW diagnostics CK
     989          ! & PTOPLWADAERO,PSOLLWADAERO,&
     990          ! & PTOPLWAD0AERO,PSOLLWAD0AERO,&
     991          ! & PTOPLWAIAERO,PSOLLWAIAERO,&
     992          ! & PLWADAERO,& !--NL
     993          !!..end
     994          ! & ok_ade, ok_aie, ok_volcan, flag_aerosol,flag_aerosol_strat,&
     995          ! & flag_aer_feedback)
     996
     997
     998          !        print *,'RADLWSW: apres RECMWF'
     999          IF (lldebug) THEN
     1000             CALL writefield_phy('zemtd_i',ZEMTD_i,klev+1)
     1001             CALL writefield_phy('zemtu_i',ZEMTU_i,klev+1)
     1002             CALL writefield_phy('ztrso_i',ZTRSO_i,klev+1)
     1003             CALL writefield_phy('zth_i',ZTH_i,klev+1)
     1004             CALL writefield_phy('zctrso',ZCTRSO,2)
     1005             CALL writefield_phy('zcemtr',ZCEMTR,2)
     1006             CALL writefield_phy('ztrsod',ZTRSOD,1)
     1007             CALL writefield_phy('zlwfc',ZLWFC,2)
     1008             CALL writefield_phy('zlwft_i',ZLWFT_i,klev+1)
     1009             CALL writefield_phy('zswfc',ZSWFC,2)
     1010             CALL writefield_phy('zswft_i',ZSWFT_i,klev+1)
     1011             CALL writefield_phy('psfswdir',PSFSWDIR,6)
     1012             CALL writefield_phy('psfswdif',PSFSWDIF,6)
     1013             CALL writefield_phy('pfsdnn',PFSDNN,1)
     1014             CALL writefield_phy('pfsdnv',PFSDNV,1)
     1015             CALL writefield_phy('ppiza_dst',PPIZA_TOT,klev)
     1016             CALL writefield_phy('pcga_dst',PCGA_TOT,klev)
     1017             CALL writefield_phy('ptaurel_dst',PTAU_TOT,klev)
     1018             CALL writefield_phy('zflux_i',ZFLUX_i,klev+1)
     1019             CALL writefield_phy('zfluc_i',ZFLUC_i,klev+1)
     1020             CALL writefield_phy('zfsdwn_i',ZFSDWN_i,klev+1)
     1021             CALL writefield_phy('zfsup_i',ZFSUP_i,klev+1)
     1022             CALL writefield_phy('zfcdwn_i',ZFCDWN_i,klev+1)
     1023             CALL writefield_phy('zfcup_i',ZFCUP_i,klev+1)
     1024          ENDIF
     1025
     1026          ! ---------
     1027          ! ---------
     1028          ! On retablit l'ordre des niveaux lmd pour les tableaux de sortie
     1029          ! D autre part, on multiplie les resultats SW par fract pour etre coherent
     1030          ! avec l ancien rayonnement AR4. Si nuit, fract=0 donc pas de
     1031          ! rayonnement SW. (MPL 260609)
     1032          DO k=0,klev
     1033             DO i=1,klon
     1034                ZEMTD(i,k+1)  = ZEMTD_i(i,k+1)
     1035                ZEMTU(i,k+1)  = ZEMTU_i(i,k+1)
     1036                ZTRSO(i,k+1)  = ZTRSO_i(i,k+1)
     1037                ZTH(i,k+1)    = ZTH_i(i,k+1)
     1038                !        ZLWFT(i,k+1)  = ZLWFT_i(i,klev+1-k)
     1039                !        ZSWFT(i,k+1)  = ZSWFT_i(i,klev+1-k)
     1040                ZFLUP(i,k+1)  = ZFLUX_i(i,1,k+1)
     1041                ZFLDN(i,k+1)  = ZFLUX_i(i,2,k+1)
     1042                ZFLUP0(i,k+1) = ZFLUC_i(i,1,k+1)
     1043                ZFLDN0(i,k+1) = ZFLUC_i(i,2,k+1)
     1044                ZFSDN(i,k+1)  = ZFSDWN_i(i,k+1)*fract(i)
     1045                ZFSDN0(i,k+1) = ZFCDWN_i(i,k+1)*fract(i)
     1046                ZFSDNC0(i,k+1)= ZFCCDWN_i(i,k+1)*fract(i)
     1047                ZFSUP (i,k+1) = ZFSUP_i(i,k+1)*fract(i)
     1048                ZFSUP0(i,k+1) = ZFCUP_i(i,k+1)*fract(i)
     1049                ZFSUPC0(i,k+1)= ZFCCUP_i(i,k+1)*fract(i)
     1050                ZFLDNC0(i,k+1)= ZFLCCDWN_i(i,k+1)
     1051                ZFLUPC0(i,k+1)= ZFLCCUP_i(i,k+1)
     1052                IF (ok_volcan) THEN
     1053                   ZSWADAERO(i,k+1)=ZSWADAERO(i,k+1)*fract(i) !--NL
     1054                ENDIF
     1055
     1056                !   Nouveau calcul car visiblement ZSWFT et ZSWFC sont nuls dans RRTM cy32
     1057                !   en sortie de radlsw.F90 - MPL 7.01.09
     1058                ZSWFT(i,k+1)  = (ZFSDWN_i(i,k+1)-ZFSUP_i(i,k+1))*fract(i)
     1059                ZSWFT0_i(i,k+1) = (ZFCDWN_i(i,k+1)-ZFCUP_i(i,k+1))*fract(i)
     1060                !        WRITE(*,'("FSDN FSUP FCDN FCUP: ",4E12.5)') ZFSDWN_i(i,k+1),&
     1061                !        ZFSUP_i(i,k+1),ZFCDWN_i(i,k+1),ZFCUP_i(i,k+1)
     1062                ZLWFT(i,k+1) =-ZFLUX_i(i,2,k+1)-ZFLUX_i(i,1,k+1)
     1063                ZLWFT0_i(i,k+1)=-ZFLUC_i(i,2,k+1)-ZFLUC_i(i,1,k+1)
     1064                !        print *,'FLUX2 FLUX1 FLUC2 FLUC1',ZFLUX_i(i,2,k+1),&
     1065                !    & ZFLUX_i(i,1,k+1),ZFLUC_i(i,2,k+1),ZFLUC_i(i,1,k+1)
     1066             ENDDO
     1067          ENDDO
     1068
     1069          !--ajout OB
     1070          ZTOPSWADAERO(:) =ZTOPSWADAERO(:) *fract(:)
     1071          ZSOLSWADAERO(:) =ZSOLSWADAERO(:) *fract(:)
     1072          ZTOPSWAD0AERO(:)=ZTOPSWAD0AERO(:)*fract(:)
     1073          ZSOLSWAD0AERO(:)=ZSOLSWAD0AERO(:)*fract(:)
     1074          ZTOPSWAIAERO(:) =ZTOPSWAIAERO(:) *fract(:)
     1075          ZSOLSWAIAERO(:) =ZSOLSWAIAERO(:) *fract(:)
     1076          ZTOPSWCF_AERO(:,1)=ZTOPSWCF_AERO(:,1)*fract(:)
     1077          ZTOPSWCF_AERO(:,2)=ZTOPSWCF_AERO(:,2)*fract(:)
     1078          ZTOPSWCF_AERO(:,3)=ZTOPSWCF_AERO(:,3)*fract(:)
     1079          ZSOLSWCF_AERO(:,1)=ZSOLSWCF_AERO(:,1)*fract(:)
     1080          ZSOLSWCF_AERO(:,2)=ZSOLSWCF_AERO(:,2)*fract(:)
     1081          ZSOLSWCF_AERO(:,3)=ZSOLSWCF_AERO(:,3)*fract(:)
     1082
     1083          ! ---------
     1084          ! ---------
     1085          ! On renseigne les champs LMDz, pour avoir la meme chose qu'en sortie de
     1086          ! LW_LMDAR4 et SW_LMDAR4
     1087
     1088          !--fraction of diffuse radiation in surface SW downward radiation
     1089          DO i = 1, kdlon
     1090             IF (fract(i).GT.0.0) THEN
     1091                zdir=SUM(PSFSWDIR(i,:))
     1092                zdif=SUM(PSFSWDIF(i,:))
     1093                zsolswfdiff(i) = zdif/(zdir+zdif)
     1094             ELSE  !--night
     1095                zsolswfdiff(i) = 1.0
     1096             ENDIF
     1097          ENDDO
     1098          !
     1099          DO i = 1, kdlon
     1100             zsolsw(i)    = ZSWFT(i,1)
     1101             zsolsw0(i)   = ZSWFT0_i(i,1)
     1102             !        zsolsw0(i)   = ZFSDN0(i,1)     -ZFSUP0(i,1)
     1103             ztopsw(i)    = ZSWFT(i,klev+1)
     1104             ztopsw0(i)   = ZSWFT0_i(i,klev+1)
     1105             !        ztopsw0(i)   = ZFSDN0(i,klev+1)-ZFSUP0(i,klev+1)
     1106             !         
     1107             !        zsollw(i)    = ZFLDN(i,1)      -ZFLUP(i,1)
     1108             !        zsollw0(i)   = ZFLDN0(i,1)     -ZFLUP0(i,1)
     1109             !        ztoplw(i)    = ZFLDN(i,klev+1) -ZFLUP(i,klev+1)
     1110             !        ztoplw0(i)   = ZFLDN0(i,klev+1)-ZFLUP0(i,klev+1)
     1111             zsollw(i)    = ZLWFT(i,1)
     1112             zsollw0(i)   = ZLWFT0_i(i,1)
     1113             ztoplw(i)    = ZLWFT(i,klev+1)*(-1)
     1114             ztoplw0(i)   = ZLWFT0_i(i,klev+1)*(-1)
     1115             !         
     1116             IF (fract(i) == 0.) THEN
    11171117!!!!! A REVOIR MPL (20090630) ca n a pas de sens quand fract=0
    1118 ! pas plus que dans le sw_AR4
    1119           zalbpla(i)   = 1.0e+39
    1120          ELSE
    1121           zalbpla(i)   = ZFSUP(i,klev+1)/ZFSDN(i,klev+1)
    1122          ENDIF
     1118                ! pas plus que dans le sw_AR4
     1119                zalbpla(i)   = 1.0e+39
     1120             ELSE
     1121                zalbpla(i)   = ZFSUP(i,klev+1)/ZFSDN(i,klev+1)
     1122             ENDIF
    11231123!!! 5 juin 2015
    11241124!!! Correction MP bug RRTM
    1125          zsollwdown(i)= -1.*ZFLDN(i,1)
    1126       ENDDO
    1127 !     print*,'OK2'
    1128 
    1129 !--add VOLMIP (surf cool or strat heat activate)
    1130       IF (flag_volc_surfstrat > 0) THEN
    1131          DO i = 1, kdlon
    1132             zsolsw(i)    = volmip_solsw(i)*fract(i)
    1133          ENDDO
    1134       ENDIF
    1135 
    1136 ! extrait de SW_AR4
    1137 !     DO k = 1, KFLEV
    1138 !        kpl1 = k+1
    1139 !        DO i = 1, KDLON
    1140 !           PHEAT(i,k) = -(ZFSUP(i,kpl1)-ZFSUP(i,k)) -(ZFSDN(i,k)-ZFSDN(i,kpl1))
    1141 !           PHEAT(i,k) = PHEAT(i,k) * RDAY*RG/RCPD / PDP(i,k)
    1142 ! ZLWFT(klon,k),ZSWFT
    1143 
    1144       DO k=1,kflev
    1145          DO i=1,kdlon
    1146            zheat(i,k)=(ZSWFT(i,k+1)-ZSWFT(i,k))*RDAY*RG/RCPD/PDP(i,k)
    1147            zheat0(i,k)=(ZSWFT0_i(i,k+1)-ZSWFT0_i(i,k))*RDAY*RG/RCPD/PDP(i,k)
    1148            zcool(i,k)=(ZLWFT(i,k)-ZLWFT(i,k+1))*RDAY*RG/RCPD/PDP(i,k)
    1149            zcool0(i,k)=(ZLWFT0_i(i,k)-ZLWFT0_i(i,k+1))*RDAY*RG/RCPD/PDP(i,k)
    1150            IF (ok_volcan) THEN
    1151               zheat_volc(i,k)=(ZSWADAERO(i,k+1)-ZSWADAERO(i,k))*RG/RCPD/PDP(i,k) !NL
    1152               zcool_volc(i,k)=(ZLWADAERO(i,k)-ZLWADAERO(i,k+1))*RG/RCPD/PDP(i,k) !NL
    1153            ENDIF
    1154 !          print *,'heat cool heat0 cool0 ',zheat(i,k),zcool(i,k),zheat0(i,k),zcool0(i,k)
    1155 !          ZFLUCUP_i(i,k)=ZFLUC_i(i,1,k)
    1156 !          ZFLUCDWN_i(i,k)=ZFLUC_i(i,2,k)         
    1157          ENDDO
    1158       ENDDO
     1125             zsollwdown(i)= -1.*ZFLDN(i,1)
     1126          ENDDO
     1127          !     print*,'OK2'
     1128
     1129          !--add VOLMIP (surf cool or strat heat activate)
     1130          IF (flag_volc_surfstrat > 0) THEN
     1131             DO i = 1, kdlon
     1132                zsolsw(i)    = volmip_solsw(i)*fract(i)
     1133             ENDDO
     1134          ENDIF
     1135
     1136          ! extrait de SW_AR4
     1137          !     DO k = 1, KFLEV
     1138          !        kpl1 = k+1
     1139          !        DO i = 1, KDLON
     1140          !           PHEAT(i,k) = -(ZFSUP(i,kpl1)-ZFSUP(i,k)) -(ZFSDN(i,k)-ZFSDN(i,kpl1))
     1141          !           PHEAT(i,k) = PHEAT(i,k) * RDAY*RG/RCPD / PDP(i,k)
     1142          ! ZLWFT(klon,k),ZSWFT
     1143
     1144          DO k=1,kflev
     1145             DO i=1,kdlon
     1146                zheat(i,k)=(ZSWFT(i,k+1)-ZSWFT(i,k))*RDAY*RG/RCPD/PDP(i,k)
     1147                zheat0(i,k)=(ZSWFT0_i(i,k+1)-ZSWFT0_i(i,k))*RDAY*RG/RCPD/PDP(i,k)
     1148                zcool(i,k)=(ZLWFT(i,k)-ZLWFT(i,k+1))*RDAY*RG/RCPD/PDP(i,k)
     1149                zcool0(i,k)=(ZLWFT0_i(i,k)-ZLWFT0_i(i,k+1))*RDAY*RG/RCPD/PDP(i,k)
     1150                IF (ok_volcan) THEN
     1151                   zheat_volc(i,k)=(ZSWADAERO(i,k+1)-ZSWADAERO(i,k))*RG/RCPD/PDP(i,k) !NL
     1152                   zcool_volc(i,k)=(ZLWADAERO(i,k)-ZLWADAERO(i,k+1))*RG/RCPD/PDP(i,k) !NL
     1153                ENDIF
     1154                !          print *,'heat cool heat0 cool0 ',zheat(i,k),zcool(i,k),zheat0(i,k),zcool0(i,k)
     1155                !          ZFLUCUP_i(i,k)=ZFLUC_i(i,1,k)
     1156                !          ZFLUCDWN_i(i,k)=ZFLUC_i(i,2,k)         
     1157             ENDDO
     1158          ENDDO
    11591159#else
    1160     abort_message="You should compile with -rrtm if running with iflag_rrtm=1"
    1161     call abort_physic(modname, abort_message, 1)
     1160          abort_message="You should compile with -rrtm if running with iflag_rrtm=1"
     1161          call abort_physic(modname, abort_message, 1)
    11621162#endif
    1163 !======================================================================
    1164 ! AI fev 2021
    1165     ELSE IF(iflag_rrtm == 2) THEN
    1166     print*,'Traitement cas iflag_rrtm = ',iflag_rrtm
    1167 !    print*,'Mise a zero des flux '
     1163          !======================================================================
     1164          ! AI fev 2021
     1165       ELSE IF(iflag_rrtm == 2) THEN
     1166          print*,'Traitement cas iflag_rrtm = ',iflag_rrtm
     1167          !    print*,'Mise a zero des flux '
    11681168#ifdef CPP_ECRAD
    1169       DO k = 1, kflev+1
    1170       DO i = 1, kdlon
    1171         ZEMTD_i(i,k)=0.
    1172         ZEMTU_i(i,k)=0.
    1173         ZTRSO_i(i,k)=0.
    1174         ZTH_i(i,k)=0.
    1175         ZLWFT_i(i,k)=0.
    1176         ZSWFT_i(i,k)=0.
    1177         ZFLUX_i(i,1,k)=0.
    1178         ZFLUX_i(i,2,k)=0.
    1179         ZFLUC_i(i,1,k)=0.
    1180         ZFLUC_i(i,2,k)=0.
    1181         ZFSDWN_i(i,k)=0.
    1182         ZFCDWN_i(i,k)=0.
    1183         ZFCCDWN_i(i,k)=0.
    1184         ZFSUP_i(i,k)=0.
    1185         ZFCUP_i(i,k)=0.
    1186         ZFCCUP_i(i,k)=0.
    1187         ZFLCCDWN_i(i,k)=0.
    1188         ZFLCCUP_i(i,k)=0.
    1189       ENDDO
    1190       ENDDO
    1191 !
    1192 ! AI ATTENTION Aerosols A REVOIR
    1193       DO i = 1, kdlon
    1194       DO k = 1, kflev
    1195       DO kk= 1, naero_spc
    1196 !      DO kk=1, NSW
    1197 !
    1198 !      PTAU_TOT(i,kflev+1-k,kk)=tau_aero_sw_rrtm(i,k,2,kk)
    1199 !      PPIZA_TOT(i,kflev+1-k,kk)=piz_aero_sw_rrtm(i,k,2,kk)
    1200 !      PCGA_TOT(i,kflev+1-k,kk)=cg_aero_sw_rrtm(i,k,2,kk)
    1201 !
    1202 !      PTAU_NAT(i,kflev+1-k,kk)=tau_aero_sw_rrtm(i,k,1,kk)
    1203 !      PPIZA_NAT(i,kflev+1-k,kk)=piz_aero_sw_rrtm(i,k,1,kk)
    1204 !      PCGA_NAT(i,kflev+1-k,kk)=cg_aero_sw_rrtm(i,k,1,kk)
    1205 !       ZAEROSOL(i,kflev+1-k,kk)=m_allaer(i,k,kk)
    1206        ZAEROSOL(i,kflev+1-k,kk)=m_allaer(i,k,kk)
    1207 !
    1208       ENDDO
    1209       ENDDO
    1210       ENDDO
    1211 !-end OB
    1212 !
    1213 !      DO i = 1, kdlon
    1214 !      DO k = 1, kflev
    1215 !      DO kk=1, NLW
    1216 !
    1217 !      PTAU_LW_TOT(i,kflev+1-k,kk)=tau_aero_lw_rrtm(i,k,2,kk)
    1218 !      PTAU_LW_NAT(i,kflev+1-k,kk)=tau_aero_lw_rrtm(i,k,1,kk)
    1219 !
    1220 !      ENDDO
    1221 !      ENDDO
    1222 !      ENDDO
    1223 !-end C. Kleinschmitt
    1224 !     
    1225       DO i = 1, kdlon
    1226       ZCTRSO(i,1)=0.
    1227       ZCTRSO(i,2)=0.
    1228       ZCEMTR(i,1)=0.
    1229       ZCEMTR(i,2)=0.
    1230       ZTRSOD(i)=0.
    1231       ZLWFC(i,1)=0.
    1232       ZLWFC(i,2)=0.
    1233       ZSWFC(i,1)=0.
    1234       ZSWFC(i,2)=0.
    1235       PFSDNN(i)=0.
    1236       PFSDNV(i)=0.
    1237       DO kk = 1, NSW
    1238         PSFSWDIR(i,kk)=0.
    1239         PSFSWDIF(i,kk)=0.
    1240       ENDDO
    1241       ENDDO
    1242 !----- Fin des mises a zero des tableaux output -------------------             
    1243 
    1244 ! On met les donnees dans l'ordre des niveaux ecrad
    1245 !         print*,'On inverse sur la verticale '
    1246          paprs_i(:,1)=paprs(:,klev+1)
    1247          DO k=1,klev
    1248             paprs_i(1:klon,k+1) =paprs(1:klon,klev+1-k)
    1249             pplay_i(1:klon,k)   =pplay(1:klon,klev+1-k)
    1250             cldfra_i(1:klon,k)  =cldfra(1:klon,klev+1-k)
    1251             PDP_i(1:klon,k)     =PDP(1:klon,klev+1-k)
    1252             t_i(1:klon,k)       =t(1:klon,klev+1-k)
    1253             q_i(1:klon,k)       =q(1:klon,klev+1-k)
    1254             qsat_i(1:klon,k)    =qsat(1:klon,klev+1-k)
    1255             flwc_i(1:klon,k)    =flwc(1:klon,klev+1-k)
    1256             fiwc_i(1:klon,k)    =fiwc(1:klon,klev+1-k)
    1257             ref_liq_i(1:klon,k) =ref_liq(1:klon,klev+1-k)*1.0e-6
    1258             ref_ice_i(1:klon,k) =ref_ice(1:klon,klev+1-k)*1.0e-6
    1259 !-OB
    1260             ref_liq_pi_i(1:klon,k) =ref_liq_pi(1:klon,klev+1-k)
    1261             ref_ice_pi_i(1:klon,k) =ref_ice_pi(1:klon,klev+1-k)
    1262          ENDDO
    1263          DO k=1,kflev
    1264             POZON_i(1:klon,k,:)=POZON(1:klon,kflev+1-k,:)
    1265 !            ZO3_DP_i(1:klon,k)=ZO3_DP(1:klon,kflev+1-k)
    1266 !            DO i=1,6
    1267             PAER_i(1:klon,k,:)=PAER(1:klon,kflev+1-k,:)
    1268 !            ENDDO
    1269          ENDDO
    1270 
    1271 ! AI 11.2021
    1272 ! Calcul de ZTH_i (temp aux interfaces 1:klev+1)
    1273 ! IFS currently sets the half-level temperature at the surface to be
    1274 ! equal to the skin temperature. The radiation scheme takes as input
    1275 ! only the half-level temperatures and assumes the Planck function to
    1276 ! vary linearly in optical depth between half levels. In the lowest
    1277 ! atmospheric layer, where the atmospheric temperature can be much
    1278 ! cooler than the skin temperature, this can lead to significant
    1279 ! differences between the effective temperature of this lowest layer
    1280 ! and the true value in the model.
    1281 ! We may approximate the temperature profile in the lowest model level
    1282 ! as piecewise linear between the top of the layer T[k-1/2], the
    1283 ! centre of the layer T[k] and the base of the layer Tskin.  The mean
    1284 ! temperature of the layer is then 0.25*T[k-1/2] + 0.5*T[k] +
    1285 ! 0.25*Tskin, which can be achieved by setting the atmospheric
    1286 ! temperature at the half-level corresponding to the surface as
    1287 ! follows:
    1288 ! AI ATTENTION fais dans interface radlw
    1289 !thermodynamics%temperature_hl(KIDIA:KFDIA,KLEV+1) &
    1290 !     &  = PTEMPERATURE(KIDIA:KFDIA,KLEV) &
    1291 !     &  + 0.5_JPRB * (PTEMPERATURE_H(KIDIA:KFDIA,KLEV+1) &
    1292 !     &               -PTEMPERATURE_H(KIDIA:KFDIA,KLEV))
    1293 
    1294          DO K=2,KLEV
    1295           DO i = 1, kdlon
    1296             ZTH_i(i,K)=&
    1297               & (t_i(i,K-1)*pplay_i(i,K-1)*(pplay_i(i,K)-paprs_i(i,K))&
    1298               & +t_i(i,K)*pplay_i(i,K)*(paprs_i(i,K)-pplay_i(i,K-1)))&
    1299               & *(1.0/(paprs_i(i,K)*(pplay_i(i,K)-pplay_i(i,K-1))))
    1300            ENDDO
    1301          ENDDO
    1302          DO i = 1, kdlon
    1303 ! Sommet
    1304             ZTH_i(i,1)=t_i(i,1)-pplay_i(i,1)*(t_i(i,1)-ZTH_i(i,2))&
    1305                       & /(pplay_i(i,1)-paprs_i(i,2))
    1306 ! Vers le sol
    1307             ZTH_i(i,KLEV+1)=t_i(i,KLEV) + 0.5 * &
    1308                             (tsol(i) - ZTH_i(i,KLEV))
    1309          ENDDO
    1310 
    1311 
    1312       print *,'RADLWSW: avant RADIATION_SCHEME '
    1313    
    1314 ! AI mars 2022
    1315     SOLARIRAD = solaire/zdist/zdist
    1316 !! diagnos pour la comparaison a la version offline
     1169          DO k = 1, kflev+1
     1170             DO i = 1, kdlon
     1171                ZEMTD_i(i,k)=0.
     1172                ZEMTU_i(i,k)=0.
     1173                ZTRSO_i(i,k)=0.
     1174                ZTH_i(i,k)=0.
     1175                ZLWFT_i(i,k)=0.
     1176                ZSWFT_i(i,k)=0.
     1177                ZFLUX_i(i,1,k)=0.
     1178                ZFLUX_i(i,2,k)=0.
     1179                ZFLUC_i(i,1,k)=0.
     1180                ZFLUC_i(i,2,k)=0.
     1181                ZFSDWN_i(i,k)=0.
     1182                ZFCDWN_i(i,k)=0.
     1183                ZFCCDWN_i(i,k)=0.
     1184                ZFSUP_i(i,k)=0.
     1185                ZFCUP_i(i,k)=0.
     1186                ZFCCUP_i(i,k)=0.
     1187                ZFLCCDWN_i(i,k)=0.
     1188                ZFLCCUP_i(i,k)=0.
     1189             ENDDO
     1190          ENDDO
     1191          !
     1192          ! AI ATTENTION Aerosols A REVOIR
     1193          DO i = 1, kdlon
     1194             DO k = 1, kflev
     1195                DO kk= 1, naero_spc
     1196                   !      DO kk=1, NSW
     1197                   !
     1198                   !      PTAU_TOT(i,kflev+1-k,kk)=tau_aero_sw_rrtm(i,k,2,kk)
     1199                   !      PPIZA_TOT(i,kflev+1-k,kk)=piz_aero_sw_rrtm(i,k,2,kk)
     1200                   !      PCGA_TOT(i,kflev+1-k,kk)=cg_aero_sw_rrtm(i,k,2,kk)
     1201                   !
     1202                   !      PTAU_NAT(i,kflev+1-k,kk)=tau_aero_sw_rrtm(i,k,1,kk)
     1203                   !      PPIZA_NAT(i,kflev+1-k,kk)=piz_aero_sw_rrtm(i,k,1,kk)
     1204                   !      PCGA_NAT(i,kflev+1-k,kk)=cg_aero_sw_rrtm(i,k,1,kk)
     1205                   !       ZAEROSOL(i,kflev+1-k,kk)=m_allaer(i,k,kk)
     1206                   ZAEROSOL(i,kflev+1-k,kk)=m_allaer(i,k,kk)
     1207                   !
     1208                ENDDO
     1209             ENDDO
     1210          ENDDO
     1211          !-end OB
     1212          !
     1213          !      DO i = 1, kdlon
     1214          !      DO k = 1, kflev
     1215          !      DO kk=1, NLW
     1216          !
     1217          !      PTAU_LW_TOT(i,kflev+1-k,kk)=tau_aero_lw_rrtm(i,k,2,kk)
     1218          !      PTAU_LW_NAT(i,kflev+1-k,kk)=tau_aero_lw_rrtm(i,k,1,kk)
     1219          !
     1220          !      ENDDO
     1221          !      ENDDO
     1222          !      ENDDO
     1223          !-end C. Kleinschmitt
     1224          !     
     1225          DO i = 1, kdlon
     1226             ZCTRSO(i,1)=0.
     1227             ZCTRSO(i,2)=0.
     1228             ZCEMTR(i,1)=0.
     1229             ZCEMTR(i,2)=0.
     1230             ZTRSOD(i)=0.
     1231             ZLWFC(i,1)=0.
     1232             ZLWFC(i,2)=0.
     1233             ZSWFC(i,1)=0.
     1234             ZSWFC(i,2)=0.
     1235             PFSDNN(i)=0.
     1236             PFSDNV(i)=0.
     1237             DO kk = 1, NSW
     1238                PSFSWDIR(i,kk)=0.
     1239                PSFSWDIF(i,kk)=0.
     1240             ENDDO
     1241          ENDDO
     1242          !----- Fin des mises a zero des tableaux output -------------------             
     1243
     1244          ! On met les donnees dans l'ordre des niveaux ecrad
     1245          !         print*,'On inverse sur la verticale '
     1246          paprs_i(:,1)=paprs(:,klev+1)
     1247          DO k=1,klev
     1248             paprs_i(1:klon,k+1) =paprs(1:klon,klev+1-k)
     1249             pplay_i(1:klon,k)   =pplay(1:klon,klev+1-k)
     1250             cldfra_i(1:klon,k)  =cldfra(1:klon,klev+1-k)
     1251             PDP_i(1:klon,k)     =PDP(1:klon,klev+1-k)
     1252             t_i(1:klon,k)       =t(1:klon,klev+1-k)
     1253             q_i(1:klon,k)       =q(1:klon,klev+1-k)
     1254             qsat_i(1:klon,k)    =qsat(1:klon,klev+1-k)
     1255             flwc_i(1:klon,k)    =flwc(1:klon,klev+1-k)
     1256             fiwc_i(1:klon,k)    =fiwc(1:klon,klev+1-k)
     1257             ref_liq_i(1:klon,k) =ref_liq(1:klon,klev+1-k)*1.0e-6
     1258             ref_ice_i(1:klon,k) =ref_ice(1:klon,klev+1-k)*1.0e-6
     1259             !-OB
     1260             ref_liq_pi_i(1:klon,k) =ref_liq_pi(1:klon,klev+1-k)
     1261             ref_ice_pi_i(1:klon,k) =ref_ice_pi(1:klon,klev+1-k)
     1262          ENDDO
     1263          DO k=1,kflev
     1264             POZON_i(1:klon,k,:)=POZON(1:klon,kflev+1-k,:)
     1265             !            ZO3_DP_i(1:klon,k)=ZO3_DP(1:klon,kflev+1-k)
     1266             !            DO i=1,6
     1267             PAER_i(1:klon,k,:)=PAER(1:klon,kflev+1-k,:)
     1268             !            ENDDO
     1269          ENDDO
     1270
     1271          ! AI 11.2021
     1272          ! Calcul de ZTH_i (temp aux interfaces 1:klev+1)
     1273          ! IFS currently sets the half-level temperature at the surface to be
     1274          ! equal to the skin temperature. The radiation scheme takes as input
     1275          ! only the half-level temperatures and assumes the Planck function to
     1276          ! vary linearly in optical depth between half levels. In the lowest
     1277          ! atmospheric layer, where the atmospheric temperature can be much
     1278          ! cooler than the skin temperature, this can lead to significant
     1279          ! differences between the effective temperature of this lowest layer
     1280          ! and the true value in the model.
     1281          ! We may approximate the temperature profile in the lowest model level
     1282          ! as piecewise linear between the top of the layer T[k-1/2], the
     1283          ! centre of the layer T[k] and the base of the layer Tskin.  The mean
     1284          ! temperature of the layer is then 0.25*T[k-1/2] + 0.5*T[k] +
     1285          ! 0.25*Tskin, which can be achieved by setting the atmospheric
     1286          ! temperature at the half-level corresponding to the surface as
     1287          ! follows:
     1288          ! AI ATTENTION fais dans interface radlw
     1289          !thermodynamics%temperature_hl(KIDIA:KFDIA,KLEV+1) &
     1290          !     &  = PTEMPERATURE(KIDIA:KFDIA,KLEV) &
     1291          !     &  + 0.5_JPRB * (PTEMPERATURE_H(KIDIA:KFDIA,KLEV+1) &
     1292          !     &               -PTEMPERATURE_H(KIDIA:KFDIA,KLEV))
     1293
     1294          DO K=2,KLEV
     1295             DO i = 1, kdlon
     1296                ZTH_i(i,K)=&
     1297                     & (t_i(i,K-1)*pplay_i(i,K-1)*(pplay_i(i,K)-paprs_i(i,K))&
     1298                     & +t_i(i,K)*pplay_i(i,K)*(paprs_i(i,K)-pplay_i(i,K-1)))&
     1299                     & *(1.0/(paprs_i(i,K)*(pplay_i(i,K)-pplay_i(i,K-1))))
     1300             ENDDO
     1301          ENDDO
     1302          DO i = 1, kdlon
     1303             ! Sommet
     1304             ZTH_i(i,1)=t_i(i,1)-pplay_i(i,1)*(t_i(i,1)-ZTH_i(i,2))&
     1305                  & /(pplay_i(i,1)-paprs_i(i,2))
     1306             ! Vers le sol
     1307             ZTH_i(i,KLEV+1)=t_i(i,KLEV) + 0.5 * &
     1308                  (tsol(i) - ZTH_i(i,KLEV))
     1309          ENDDO
     1310
     1311
     1312          print *,'RADLWSW: avant RADIATION_SCHEME '
     1313
     1314          ! AI mars 2022
     1315          SOLARIRAD = solaire/zdist/zdist
     1316          !! diagnos pour la comparaison a la version offline
    13171317!!! - Gas en VMR pour offline et MMR pour online
    13181318!!! - on utilise pour solarirrad une valeur constante
    1319     if (lldebug_for_offline) then
    1320        SOLARIRAD = 1366.0896
    1321        ZCH4_off = CH4_ppb*1e-9
    1322        ZN2O_off = N2O_ppb*1e-9
    1323        ZNO2_off = 0.0
    1324        ZCFC11_off = CFC11_ppt*1e-12
    1325        ZCFC12_off = CFC12_ppt*1e-12
    1326        ZHCFC22_off = 0.0
    1327        ZCCL4_off = 0.0
    1328        ZO2_off = 0.0
    1329        ZCO2_off = co2_ppm*1e-6
    1330 
    1331         CALL writefield_phy('rmu0',rmu0,1)
    1332         CALL writefield_phy('tsol',tsol,1)
    1333         CALL writefield_phy('emissiv_out',ZEMIS,1)
    1334         CALL writefield_phy('paprs_i',paprs_i,klev+1)
    1335         CALL writefield_phy('ZTH_i',ZTH_i,klev+1)
    1336         CALL writefield_phy('cldfra_i',cldfra_i,klev)
    1337         CALL writefield_phy('q_i',q_i,klev)
    1338         CALL writefield_phy('fiwc_i',fiwc_i,klev)
    1339         CALL writefield_phy('flwc_i',flwc_i,klev)
    1340         CALL writefield_phy('palbd_new',PALBD_NEW,NSW)
    1341         CALL writefield_phy('palbp_new',PALBP_NEW,NSW)
    1342         CALL writefield_phy('POZON',POZON_i(:,:,1),klev)
    1343         CALL writefield_phy('ZCO2',ZCO2_off,klev)
    1344         CALL writefield_phy('ZCH4',ZCH4_off,klev)
    1345         CALL writefield_phy('ZN2O',ZN2O_off,klev)
    1346         CALL writefield_phy('ZO2',ZO2_off,klev)
    1347         CALL writefield_phy('ZNO2',ZNO2_off,klev)
    1348         CALL writefield_phy('ZCFC11',ZCFC11_off,klev)
    1349         CALL writefield_phy('ZCFC12',ZCFC12_off,klev)
    1350         CALL writefield_phy('ZHCFC22',ZHCFC22_off,klev)
    1351         CALL writefield_phy('ZCCL4',ZCCL4_off,klev)
    1352         CALL writefield_phy('ref_liq_i',ref_liq_i,klev)
    1353         CALL writefield_phy('ref_ice_i',ref_ice_i,klev)
    1354       endif
    1355 ! lldebug_for_offline
    1356 
    1357   if (namelist_ecrad_file.eq.'namelist_ecrad') then
    1358       print*,' 1er apell Ecrad : ok_3Deffect, namelist_ecrad_file = ', &
    1359           ok_3Deffect, namelist_ecrad_file   
    1360       CALL RADIATION_SCHEME &
    1361       & (ist, iend, klon, klev, naero_spc, NSW, &
    1362       & namelist_ecrad_file, ok_3Deffect, &
    1363       & debut, ok_volcan, flag_aerosol_strat, &
    1364       & day_cur, current_time, &
    1365 !       Cste solaire/(d_Terre-Soleil)**2
    1366       & SOLARIRAD, &
    1367 !       Cos(angle zin), temp sol             
    1368       & rmu0, tsol, &
    1369 !       Albedo diffuse et directe
    1370       & PALBD_NEW,PALBP_NEW, &   
    1371 !       Emessivite : PEMIS_WINDOW (???), &
    1372       & ZEMIS, ZEMISW, &
    1373 !       longitude(rad), sin(latitude), PMASQ_ ???
    1374       & ZGELAM, ZGEMU, &
    1375 !       Temp et pres aux interf, vapeur eau, Satur spec humid
    1376       & paprs_i, ZTH_i, q_i, qsat_i, &
    1377 !       Gas
    1378        & ZCO2, ZCH4, ZN2O, ZNO2, ZCFC11, ZCFC12, ZHCFC22, &
    1379        & ZCCL4, POZON_i(:,:,1), ZO2, &
    1380 !       nuages :
    1381       & cldfra_i, flwc_i, fiwc_i, ZQ_SNOW, &
    1382 !       rayons effectifs des gouttelettes             
    1383       & ref_liq_i, ref_ice_i, &
    1384 !       aerosols
    1385       & ZAEROSOL_OLD, ZAEROSOL, &
    1386 ! Outputs
    1387 !       Net flux :
    1388       & ZSWFT_i, ZLWFT_i, ZSWFT0_ii, ZLWFT0_ii, &
    1389 !       DWN flux :
    1390       & ZFSDWN_i, ZFLUX_i(:,2,:), ZFCDWN_i, ZFLUC_i(:,2,:), &
    1391 !       UP flux :
    1392       & ZFSUP_i, ZFLUX_i(:,1,:), ZFCUP_i, ZFLUC_i(:,1,:), &
    1393 !       Surf Direct flux : ATTENTION
    1394       & ZFLUX_DIR, ZFLUX_DIR_CLEAR, ZFLUX_DIR_INTO_SUN, &
    1395 !       UV and para flux
    1396       & ZFLUX_UV, ZFLUX_PAR, ZFLUX_PAR_CLEAR, &
    1397 !      & ZFLUX_SW_DN_TOA,
    1398       & ZEMIS_OUT, ZLWDERIVATIVE, &
    1399       & PSFSWDIF, PSFSWDIR, &
    1400       & cloud_cover_sw)
    1401   else
    1402    print*,' 2e apell Ecrad : ok_3Deffect, namelist_ecrad_file = ', &
    1403           ok_3Deffect, namelist_ecrad_file       
    1404    CALL RADIATION_SCHEME_S2 &
    1405       & (ist, iend, klon, klev, naero_grp, NSW, &
    1406       & namelist_ecrad_file, ok_3Deffect, &
    1407       & debut, ok_volcan, flag_aerosol_strat, &
    1408       & day_cur, current_time, &
    1409 !       Cste solaire/(d_Terre-Soleil)**2
    1410       & SOLARIRAD, &
    1411 !       Cos(angle zin), temp sol             
    1412       & rmu0, tsol, &
    1413 !       Albedo diffuse et directe
    1414       & PALBD_NEW,PALBP_NEW, &
    1415 !       Emessivite : PEMIS_WINDOW (???), &
    1416       & ZEMIS, ZEMISW, &
    1417 !       longitude(rad), sin(latitude), PMASQ_ ???
    1418       & ZGELAM, ZGEMU, &
    1419 !       Temp et pres aux interf, vapeur eau, Satur spec humid
    1420       & paprs_i, ZTH_i, q_i, qsat_i, &
    1421 !       Gas
    1422        & ZCO2, ZCH4, ZN2O, ZNO2, ZCFC11, ZCFC12, ZHCFC22, &
    1423        & ZCCL4, POZON_i(:,:,1), ZO2, &
    1424 !       nuages :
    1425       & cldfra_i, flwc_i, fiwc_i, ZQ_SNOW, &
    1426 !       rayons effectifs des gouttelettes             
    1427       & ref_liq_i, ref_ice_i, &
    1428 !       aerosols
    1429      & ZAEROSOL_OLD, ZAEROSOL, &
    1430 ! Outputs
    1431 !       Net flux :
    1432       & ZSWFT_i, ZLWFT_i, ZSWFT0_ii, ZLWFT0_ii, &
    1433 !       DWN flux :
    1434       & ZFSDWN_i, ZFLUX_i(:,2,:), ZFCDWN_i, ZFLUC_i(:,2,:), &
    1435 !       UP flux :
    1436       & ZFSUP_i, ZFLUX_i(:,1,:), ZFCUP_i, ZFLUC_i(:,1,:), &
    1437 !       Surf Direct flux : ATTENTION
    1438       & ZFLUX_DIR, ZFLUX_DIR_CLEAR, ZFLUX_DIR_INTO_SUN, &
    1439 !       UV and para flux
    1440       & ZFLUX_UV, ZFLUX_PAR, ZFLUX_PAR_CLEAR, &
    1441 !      & ZFLUX_SW_DN_TOA,
    1442       & ZEMIS_OUT, ZLWDERIVATIVE, &
    1443       & PSFSWDIF, PSFSWDIR, &
    1444       & cloud_cover_sw)
    1445   endif
    1446 
    1447 
    1448       print *,'========= RADLWSW: apres RADIATION_SCHEME ==================== '
    1449 
    1450      if (lldebug_for_offline) then
    1451         CALL writefield_phy('FLUX_LW',ZLWFT_i,klev+1)
    1452         CALL writefield_phy('FLUX_LW_CLEAR',ZLWFT0_ii,klev+1)
    1453         CALL writefield_phy('FLUX_SW',ZSWFT_i,klev+1)
    1454         CALL writefield_phy('FLUX_SW_CLEAR',ZSWFT0_ii,klev+1)
    1455         CALL writefield_phy('FLUX_DN_SW',ZFSDWN_i,klev+1)
    1456         CALL writefield_phy('FLUX_DN_LW',ZFLUX_i(:,2,:),klev+1)
    1457         CALL writefield_phy('FLUX_DN_SW_CLEAR',ZFCDWN_i,klev+1)
    1458         CALL writefield_phy('FLUX_DN_LW_CLEAR',ZFLUC_i(:,2,:),klev+1)
    1459         CALL writefield_phy('PSFSWDIR',PSFSWDIR,6)
    1460         CALL writefield_phy('PSFSWDIF',PSFSWDIF,6)
    1461         CALL writefield_phy('FLUX_UP_LW',ZFLUX_i(:,1,:),klev+1)
    1462         CALL writefield_phy('FLUX_UP_LW_CLEAR',ZFLUC_i(:,1,:),klev+1)
    1463         CALL writefield_phy('FLUX_UP_SW',ZFSUP_i,klev+1)
    1464         CALL writefield_phy('FLUX_UP_SW_CLEAR',ZFCUP_i,klev+1)
    1465       endif
    1466 
    1467 ! ---------
    1468 ! On retablit l'ordre des niveaux lmd pour les tableaux de sortie
    1469 ! D autre part, on multiplie les resultats SW par fract pour etre coherent
    1470 ! avec l ancien rayonnement AR4. Si nuit, fract=0 donc pas de
    1471 ! rayonnement SW. (MPL 260609)
    1472       print*,'On retablit l ordre des niveaux verticaux pour LMDZ'
    1473       print*,'On multiplie les flux SW par fract et LW dwn par -1'
    1474       DO k=0,klev
    1475          DO i=1,klon
    1476          ZEMTD(i,k+1)  = ZEMTD_i(i,klev+1-k)
    1477          ZEMTU(i,k+1)  = ZEMTU_i(i,klev+1-k)
    1478          ZTRSO(i,k+1)  = ZTRSO_i(i,klev+1-k)
    1479 !         ZTH(i,k+1)    = ZTH_i(i,klev+1-k)
    1480 ! AI ATTENTION
    1481           ZLWFT(i,k+1)  = ZLWFT_i(i,klev+1-k)
    1482           ZSWFT(i,k+1)  = ZSWFT_i(i,klev+1-k)*fract(i)
    1483           ZSWFT0_i(i,k+1) = ZSWFT0_ii(i,klev+1-k)*fract(i)
    1484           ZLWFT0_i(i,k+1) = ZLWFT0_ii(i,klev+1-k)
    1485 !
    1486          ZFLUP(i,k+1)  = ZFLUX_i(i,1,klev+1-k)
    1487          ZFLDN(i,k+1)  = -1.*ZFLUX_i(i,2,klev+1-k)
    1488          ZFLUP0(i,k+1) = ZFLUC_i(i,1,klev+1-k)
    1489          ZFLDN0(i,k+1) = -1.*ZFLUC_i(i,2,klev+1-k)
    1490          ZFSDN(i,k+1)  = ZFSDWN_i(i,klev+1-k)*fract(i)
    1491          ZFSDN0(i,k+1) = ZFCDWN_i(i,klev+1-k)*fract(i)
    1492          ZFSDNC0(i,k+1)= ZFCCDWN_i(i,klev+1-k)*fract(i)
    1493          ZFSUP (i,k+1) = ZFSUP_i(i,klev+1-k)*fract(i)
    1494          ZFSUP0(i,k+1) = ZFCUP_i(i,klev+1-k)*fract(i)
    1495          ZFSUPC0(i,k+1)= ZFCCUP_i(i,klev+1-k)*fract(i)
    1496          ZFLDNC0(i,k+1)= -1.*ZFLCCDWN_i(i,klev+1-k)
    1497          ZFLUPC0(i,k+1)= ZFLCCUP_i(i,klev+1-k)
    1498          IF (ok_volcan) THEN
    1499             ZSWADAERO(i,k+1)=ZSWADAERO(i,klev+1-k)*fract(i) !--NL
    1500          ENDIF
    1501          
    1502 !   Nouveau calcul car visiblement ZSWFT et ZSWFC sont nuls dans RRTM cy32
    1503 !   en sortie de radlsw.F90 - MPL 7.01.09
    1504 ! AI ATTENTION
    1505 !         ZSWFT(i,k+1)  = (ZFSDWN_i(i,k+1)-ZFSUP_i(i,k+1))*fract(i)
    1506 !         ZSWFT0_i(i,k+1) = (ZFCDWN_i(i,k+1)-ZFCUP_i(i,k+1))*fract(i)
    1507 !         ZLWFT(i,k+1) =-ZFLUX_i(i,2,k+1)-ZFLUX_i(i,1,k+1)
    1508 !         ZLWFT0_i(i,k+1)=-ZFLUC_i(i,2,k+1)-ZFLUC_i(i,1,k+1)
    1509          ENDDO
    1510       ENDDO
    1511 
    1512 !--ajout OB
    1513       ZTOPSWADAERO(:) =ZTOPSWADAERO(:) *fract(:)
    1514       ZSOLSWADAERO(:) =ZSOLSWADAERO(:) *fract(:)
    1515       ZTOPSWAD0AERO(:)=ZTOPSWAD0AERO(:)*fract(:)
    1516       ZSOLSWAD0AERO(:)=ZSOLSWAD0AERO(:)*fract(:)
    1517       ZTOPSWAIAERO(:) =ZTOPSWAIAERO(:) *fract(:)
    1518       ZSOLSWAIAERO(:) =ZSOLSWAIAERO(:) *fract(:)
    1519       ZTOPSWCF_AERO(:,1)=ZTOPSWCF_AERO(:,1)*fract(:)
    1520       ZTOPSWCF_AERO(:,2)=ZTOPSWCF_AERO(:,2)*fract(:)
    1521       ZTOPSWCF_AERO(:,3)=ZTOPSWCF_AERO(:,3)*fract(:)
    1522       ZSOLSWCF_AERO(:,1)=ZSOLSWCF_AERO(:,1)*fract(:)
    1523       ZSOLSWCF_AERO(:,2)=ZSOLSWCF_AERO(:,2)*fract(:)
    1524       ZSOLSWCF_AERO(:,3)=ZSOLSWCF_AERO(:,3)*fract(:)
    1525 
    1526 ! ---------
    1527 ! On renseigne les champs LMDz, pour avoir la meme chose qu'en sortie de
    1528 ! LW_LMDAR4 et SW_LMDAR4
    1529 
    1530       !--fraction of diffuse radiation in surface SW downward radiation
    1531       DO i = 1, kdlon
    1532          zdir=SUM(PSFSWDIR(i,:))
    1533          zdif=SUM(PSFSWDIF(i,:))
    1534        IF (fract(i).GT.0.0.and.(zdir+zdif).gt.seuilmach) THEN
    1535          zsolswfdiff(i) = zdif/(zdir+zdif)
    1536        ELSE  !--night
    1537          zsolswfdiff(i) = 1.0
     1319          if (lldebug_for_offline) then
     1320             SOLARIRAD = 1366.0896
     1321             ZCH4_off = CH4_ppb*1e-9
     1322             ZN2O_off = N2O_ppb*1e-9
     1323             ZNO2_off = 0.0
     1324             ZCFC11_off = CFC11_ppt*1e-12
     1325             ZCFC12_off = CFC12_ppt*1e-12
     1326             ZHCFC22_off = 0.0
     1327             ZCCL4_off = 0.0
     1328             ZO2_off = 0.0
     1329             ZCO2_off = co2_ppm*1e-6
     1330
     1331             CALL writefield_phy('rmu0',rmu0,1)
     1332             CALL writefield_phy('tsol',tsol,1)
     1333             CALL writefield_phy('emissiv_out',ZEMIS,1)
     1334             CALL writefield_phy('paprs_i',paprs_i,klev+1)
     1335             CALL writefield_phy('ZTH_i',ZTH_i,klev+1)
     1336             CALL writefield_phy('cldfra_i',cldfra_i,klev)
     1337             CALL writefield_phy('q_i',q_i,klev)
     1338             CALL writefield_phy('fiwc_i',fiwc_i,klev)
     1339             CALL writefield_phy('flwc_i',flwc_i,klev)
     1340             CALL writefield_phy('palbd_new',PALBD_NEW,NSW)
     1341             CALL writefield_phy('palbp_new',PALBP_NEW,NSW)
     1342             CALL writefield_phy('POZON',POZON_i(:,:,1),klev)
     1343             CALL writefield_phy('ZCO2',ZCO2_off,klev)
     1344             CALL writefield_phy('ZCH4',ZCH4_off,klev)
     1345             CALL writefield_phy('ZN2O',ZN2O_off,klev)
     1346             CALL writefield_phy('ZO2',ZO2_off,klev)
     1347             CALL writefield_phy('ZNO2',ZNO2_off,klev)
     1348             CALL writefield_phy('ZCFC11',ZCFC11_off,klev)
     1349             CALL writefield_phy('ZCFC12',ZCFC12_off,klev)
     1350             CALL writefield_phy('ZHCFC22',ZHCFC22_off,klev)
     1351             CALL writefield_phy('ZCCL4',ZCCL4_off,klev)
     1352             CALL writefield_phy('ref_liq_i',ref_liq_i,klev)
     1353             CALL writefield_phy('ref_ice_i',ref_ice_i,klev)
     1354          endif
     1355          ! lldebug_for_offline
     1356
     1357          if (namelist_ecrad_file.eq.'namelist_ecrad') then
     1358             print*,' 1er apell Ecrad : ok_3Deffect, namelist_ecrad_file = ', &
     1359                  ok_3Deffect, namelist_ecrad_file   
     1360             CALL RADIATION_SCHEME &
     1361                  & (ist, iend, klon, klev, naero_spc, NSW, &
     1362                  & namelist_ecrad_file, ok_3Deffect, &
     1363                  & debut, ok_volcan, flag_aerosol_strat, &
     1364                  & day_cur, current_time, &
     1365                  !       Cste solaire/(d_Terre-Soleil)**2
     1366                  & SOLARIRAD, &
     1367                  !       Cos(angle zin), temp sol             
     1368                  & rmu0, tsol, &
     1369                  !       Albedo diffuse et directe
     1370                  & PALBD_NEW,PALBP_NEW, &   
     1371                  !       Emessivite : PEMIS_WINDOW (???), &
     1372                  & ZEMIS, ZEMISW, &
     1373                  !       longitude(rad), sin(latitude), PMASQ_ ???
     1374                  & ZGELAM, ZGEMU, &
     1375                  !       Temp et pres aux interf, vapeur eau, Satur spec humid
     1376                  & paprs_i, ZTH_i, q_i, qsat_i, &
     1377                  !       Gas
     1378                  & ZCO2, ZCH4, ZN2O, ZNO2, ZCFC11, ZCFC12, ZHCFC22, &
     1379                  & ZCCL4, POZON_i(:,:,1), ZO2, &
     1380                  !       nuages :
     1381                  & cldfra_i, flwc_i, fiwc_i, ZQ_SNOW, &
     1382                  !       rayons effectifs des gouttelettes             
     1383                  & ref_liq_i, ref_ice_i, &
     1384                  !       aerosols
     1385                  & ZAEROSOL_OLD, ZAEROSOL, &
     1386                  ! Outputs
     1387                  !       Net flux :
     1388                  & ZSWFT_i, ZLWFT_i, ZSWFT0_ii, ZLWFT0_ii, &
     1389                  !       DWN flux :
     1390                  & ZFSDWN_i, ZFLUX_i(:,2,:), ZFCDWN_i, ZFLUC_i(:,2,:), &
     1391                  !       UP flux :
     1392                  & ZFSUP_i, ZFLUX_i(:,1,:), ZFCUP_i, ZFLUC_i(:,1,:), &
     1393                  !       Surf Direct flux : ATTENTION
     1394                  & ZFLUX_DIR, ZFLUX_DIR_CLEAR, ZFLUX_DIR_INTO_SUN, &
     1395                  !       UV and para flux
     1396                  & ZFLUX_UV, ZFLUX_PAR, ZFLUX_PAR_CLEAR, &
     1397                  !      & ZFLUX_SW_DN_TOA,
     1398                  & ZEMIS_OUT, ZLWDERIVATIVE, &
     1399                  & PSFSWDIF, PSFSWDIR, &
     1400                  & cloud_cover_sw)
     1401          else
     1402             print*,' 2e apell Ecrad : ok_3Deffect, namelist_ecrad_file = ', &
     1403                  ok_3Deffect, namelist_ecrad_file       
     1404             CALL RADIATION_SCHEME_S2 &
     1405                  & (ist, iend, klon, klev, naero_grp, NSW, &
     1406                  & namelist_ecrad_file, ok_3Deffect, &
     1407                  & debut, ok_volcan, flag_aerosol_strat, &
     1408                  & day_cur, current_time, &
     1409                  !       Cste solaire/(d_Terre-Soleil)**2
     1410                  & SOLARIRAD, &
     1411                  !       Cos(angle zin), temp sol             
     1412                  & rmu0, tsol, &
     1413                  !       Albedo diffuse et directe
     1414                  & PALBD_NEW,PALBP_NEW, &
     1415                  !       Emessivite : PEMIS_WINDOW (???), &
     1416                  & ZEMIS, ZEMISW, &
     1417                  !       longitude(rad), sin(latitude), PMASQ_ ???
     1418                  & ZGELAM, ZGEMU, &
     1419                  !       Temp et pres aux interf, vapeur eau, Satur spec humid
     1420                  & paprs_i, ZTH_i, q_i, qsat_i, &
     1421                  !       Gas
     1422                  & ZCO2, ZCH4, ZN2O, ZNO2, ZCFC11, ZCFC12, ZHCFC22, &
     1423                  & ZCCL4, POZON_i(:,:,1), ZO2, &
     1424                  !       nuages :
     1425                  & cldfra_i, flwc_i, fiwc_i, ZQ_SNOW, &
     1426                  !       rayons effectifs des gouttelettes             
     1427                  & ref_liq_i, ref_ice_i, &
     1428                  !       aerosols
     1429                  & ZAEROSOL_OLD, ZAEROSOL, &
     1430                  ! Outputs
     1431                  !       Net flux :
     1432                  & ZSWFT_i, ZLWFT_i, ZSWFT0_ii, ZLWFT0_ii, &
     1433                  !       DWN flux :
     1434                  & ZFSDWN_i, ZFLUX_i(:,2,:), ZFCDWN_i, ZFLUC_i(:,2,:), &
     1435                  !       UP flux :
     1436                  & ZFSUP_i, ZFLUX_i(:,1,:), ZFCUP_i, ZFLUC_i(:,1,:), &
     1437                  !       Surf Direct flux : ATTENTION
     1438                  & ZFLUX_DIR, ZFLUX_DIR_CLEAR, ZFLUX_DIR_INTO_SUN, &
     1439                  !       UV and para flux
     1440                  & ZFLUX_UV, ZFLUX_PAR, ZFLUX_PAR_CLEAR, &
     1441                  !      & ZFLUX_SW_DN_TOA,
     1442                  & ZEMIS_OUT, ZLWDERIVATIVE, &
     1443                  & PSFSWDIF, PSFSWDIR, &
     1444                  & cloud_cover_sw)
     1445          endif
     1446
     1447
     1448          print *,'========= RADLWSW: apres RADIATION_SCHEME ==================== '
     1449
     1450          if (lldebug_for_offline) then
     1451             CALL writefield_phy('FLUX_LW',ZLWFT_i,klev+1)
     1452             CALL writefield_phy('FLUX_LW_CLEAR',ZLWFT0_ii,klev+1)
     1453             CALL writefield_phy('FLUX_SW',ZSWFT_i,klev+1)
     1454             CALL writefield_phy('FLUX_SW_CLEAR',ZSWFT0_ii,klev+1)
     1455             CALL writefield_phy('FLUX_DN_SW',ZFSDWN_i,klev+1)
     1456             CALL writefield_phy('FLUX_DN_LW',ZFLUX_i(:,2,:),klev+1)
     1457             CALL writefield_phy('FLUX_DN_SW_CLEAR',ZFCDWN_i,klev+1)
     1458             CALL writefield_phy('FLUX_DN_LW_CLEAR',ZFLUC_i(:,2,:),klev+1)
     1459             CALL writefield_phy('PSFSWDIR',PSFSWDIR,6)
     1460             CALL writefield_phy('PSFSWDIF',PSFSWDIF,6)
     1461             CALL writefield_phy('FLUX_UP_LW',ZFLUX_i(:,1,:),klev+1)
     1462             CALL writefield_phy('FLUX_UP_LW_CLEAR',ZFLUC_i(:,1,:),klev+1)
     1463             CALL writefield_phy('FLUX_UP_SW',ZFSUP_i,klev+1)
     1464             CALL writefield_phy('FLUX_UP_SW_CLEAR',ZFCUP_i,klev+1)
     1465          endif
     1466
     1467          ! ---------
     1468          ! On retablit l'ordre des niveaux lmd pour les tableaux de sortie
     1469          ! D autre part, on multiplie les resultats SW par fract pour etre coherent
     1470          ! avec l ancien rayonnement AR4. Si nuit, fract=0 donc pas de
     1471          ! rayonnement SW. (MPL 260609)
     1472          print*,'On retablit l ordre des niveaux verticaux pour LMDZ'
     1473          print*,'On multiplie les flux SW par fract et LW dwn par -1'
     1474          DO k=0,klev
     1475             DO i=1,klon
     1476                ZEMTD(i,k+1)  = ZEMTD_i(i,klev+1-k)
     1477                ZEMTU(i,k+1)  = ZEMTU_i(i,klev+1-k)
     1478                ZTRSO(i,k+1)  = ZTRSO_i(i,klev+1-k)
     1479                !         ZTH(i,k+1)    = ZTH_i(i,klev+1-k)
     1480                ! AI ATTENTION
     1481                ZLWFT(i,k+1)  = ZLWFT_i(i,klev+1-k)
     1482                ZSWFT(i,k+1)  = ZSWFT_i(i,klev+1-k)*fract(i)
     1483                ZSWFT0_i(i,k+1) = ZSWFT0_ii(i,klev+1-k)*fract(i)
     1484                ZLWFT0_i(i,k+1) = ZLWFT0_ii(i,klev+1-k)
     1485                !
     1486                ZFLUP(i,k+1)  = ZFLUX_i(i,1,klev+1-k)
     1487                ZFLDN(i,k+1)  = -1.*ZFLUX_i(i,2,klev+1-k)
     1488                ZFLUP0(i,k+1) = ZFLUC_i(i,1,klev+1-k)
     1489                ZFLDN0(i,k+1) = -1.*ZFLUC_i(i,2,klev+1-k)
     1490                ZFSDN(i,k+1)  = ZFSDWN_i(i,klev+1-k)*fract(i)
     1491                ZFSDN0(i,k+1) = ZFCDWN_i(i,klev+1-k)*fract(i)
     1492                ZFSDNC0(i,k+1)= ZFCCDWN_i(i,klev+1-k)*fract(i)
     1493                ZFSUP (i,k+1) = ZFSUP_i(i,klev+1-k)*fract(i)
     1494                ZFSUP0(i,k+1) = ZFCUP_i(i,klev+1-k)*fract(i)
     1495                ZFSUPC0(i,k+1)= ZFCCUP_i(i,klev+1-k)*fract(i)
     1496                ZFLDNC0(i,k+1)= -1.*ZFLCCDWN_i(i,klev+1-k)
     1497                ZFLUPC0(i,k+1)= ZFLCCUP_i(i,klev+1-k)
     1498                IF (ok_volcan) THEN
     1499                   ZSWADAERO(i,k+1)=ZSWADAERO(i,klev+1-k)*fract(i) !--NL
     1500                ENDIF
     1501
     1502                !   Nouveau calcul car visiblement ZSWFT et ZSWFC sont nuls dans RRTM cy32
     1503                !   en sortie de radlsw.F90 - MPL 7.01.09
     1504                ! AI ATTENTION
     1505                !         ZSWFT(i,k+1)  = (ZFSDWN_i(i,k+1)-ZFSUP_i(i,k+1))*fract(i)
     1506                !         ZSWFT0_i(i,k+1) = (ZFCDWN_i(i,k+1)-ZFCUP_i(i,k+1))*fract(i)
     1507                !         ZLWFT(i,k+1) =-ZFLUX_i(i,2,k+1)-ZFLUX_i(i,1,k+1)
     1508                !         ZLWFT0_i(i,k+1)=-ZFLUC_i(i,2,k+1)-ZFLUC_i(i,1,k+1)
     1509             ENDDO
     1510          ENDDO
     1511
     1512          !--ajout OB
     1513          ZTOPSWADAERO(:) =ZTOPSWADAERO(:) *fract(:)
     1514          ZSOLSWADAERO(:) =ZSOLSWADAERO(:) *fract(:)
     1515          ZTOPSWAD0AERO(:)=ZTOPSWAD0AERO(:)*fract(:)
     1516          ZSOLSWAD0AERO(:)=ZSOLSWAD0AERO(:)*fract(:)
     1517          ZTOPSWAIAERO(:) =ZTOPSWAIAERO(:) *fract(:)
     1518          ZSOLSWAIAERO(:) =ZSOLSWAIAERO(:) *fract(:)
     1519          ZTOPSWCF_AERO(:,1)=ZTOPSWCF_AERO(:,1)*fract(:)
     1520          ZTOPSWCF_AERO(:,2)=ZTOPSWCF_AERO(:,2)*fract(:)
     1521          ZTOPSWCF_AERO(:,3)=ZTOPSWCF_AERO(:,3)*fract(:)
     1522          ZSOLSWCF_AERO(:,1)=ZSOLSWCF_AERO(:,1)*fract(:)
     1523          ZSOLSWCF_AERO(:,2)=ZSOLSWCF_AERO(:,2)*fract(:)
     1524          ZSOLSWCF_AERO(:,3)=ZSOLSWCF_AERO(:,3)*fract(:)
     1525
     1526          ! ---------
     1527          ! On renseigne les champs LMDz, pour avoir la meme chose qu'en sortie de
     1528          ! LW_LMDAR4 et SW_LMDAR4
     1529
     1530          !--fraction of diffuse radiation in surface SW downward radiation
     1531          DO i = 1, kdlon
     1532             zdir=SUM(PSFSWDIR(i,:))
     1533             zdif=SUM(PSFSWDIF(i,:))
     1534             IF (fract(i).GT.0.0.and.(zdir+zdif).gt.seuilmach) THEN
     1535                zsolswfdiff(i) = zdif/(zdir+zdif)
     1536             ELSE  !--night
     1537                zsolswfdiff(i) = 1.0
     1538             ENDIF
     1539          ENDDO
     1540          !
     1541          DO i = 1, kdlon
     1542             zsolsw(i)    = ZSWFT(i,1)
     1543             zsolsw0(i)   = ZSWFT0_i(i,1)
     1544             ztopsw(i)    = ZSWFT(i,klev+1)
     1545             ztopsw0(i)   = ZSWFT0_i(i,klev+1)
     1546             zsollw(i)    = ZLWFT(i,1)
     1547             zsollw0(i)   = ZLWFT0_i(i,1)
     1548             ztoplw(i)    = ZLWFT(i,klev+1)*(-1)
     1549             ztoplw0(i)   = ZLWFT0_i(i,klev+1)*(-1)
     1550             !         
     1551             zsollwdown(i)= -1.*ZFLDN(i,1)
     1552          ENDDO
     1553
     1554          DO k=1,kflev
     1555             DO i=1,kdlon
     1556                zheat(i,k)=(ZSWFT(i,k+1)-ZSWFT(i,k))*RDAY*RG/RCPD/PDP(i,k)
     1557                zheat0(i,k)=(ZSWFT0_i(i,k+1)-ZSWFT0_i(i,k))*RDAY*RG/RCPD/PDP(i,k)
     1558                zcool(i,k)=(ZLWFT(i,k)-ZLWFT(i,k+1))*RDAY*RG/RCPD/PDP(i,k)
     1559                zcool0(i,k)=(ZLWFT0_i(i,k)-ZLWFT0_i(i,k+1))*RDAY*RG/RCPD/PDP(i,k)
     1560                IF (ok_volcan) THEN
     1561                   zheat_volc(i,k)=(ZSWADAERO(i,k+1)-ZSWADAERO(i,k))*RG/RCPD/PDP(i,k) !NL
     1562                   zcool_volc(i,k)=(ZLWADAERO(i,k)-ZLWADAERO(i,k+1))*RG/RCPD/PDP(i,k) !NL
     1563                ENDIF
     1564             ENDDO
     1565          ENDDO
     1566#endif 
     1567          print*,'Fin traitement ECRAD'
     1568          ! Fin ECRAD
     1569       ENDIF        ! iflag_rrtm
     1570       ! ecrad
     1571       !======================================================================
     1572
     1573       DO i = 1, kdlon
     1574          topsw(iof+i) = ztopsw(i)
     1575          toplw(iof+i) = ztoplw(i)
     1576          solsw(iof+i) = zsolsw(i)
     1577          solswfdiff(iof+i) = zsolswfdiff(i)
     1578          sollw(iof+i) = zsollw(i)
     1579          sollwdown(iof+i) = zsollwdown(i)
     1580          DO k = 1, kflev+1
     1581             lwdn0 ( iof+i,k)   = ZFLDN0 ( i,k)
     1582             lwdn  ( iof+i,k)   = ZFLDN  ( i,k)
     1583             lwup0 ( iof+i,k)   = ZFLUP0 ( i,k)
     1584             lwup  ( iof+i,k)   = ZFLUP  ( i,k)
     1585          ENDDO
     1586          topsw0(iof+i) = ztopsw0(i)
     1587          toplw0(iof+i) = ztoplw0(i)
     1588          solsw0(iof+i) = zsolsw0(i)
     1589          sollw0(iof+i) = zsollw0(i)
     1590          albpla(iof+i) = zalbpla(i)
     1591
     1592          DO k = 1, kflev+1
     1593             swdnc0( iof+i,k)   = ZFSDNC0( i,k)
     1594             swdn0 ( iof+i,k)   = ZFSDN0 ( i,k)
     1595             swdn  ( iof+i,k)   = ZFSDN  ( i,k)
     1596             swupc0( iof+i,k)   = ZFSUPC0( i,k)
     1597             swup0 ( iof+i,k)   = ZFSUP0 ( i,k)
     1598             swup  ( iof+i,k)   = ZFSUP  ( i,k)
     1599             lwdnc0( iof+i,k)   = ZFLDNC0( i,k)
     1600             lwupc0( iof+i,k)   = ZFLUPC0( i,k)
     1601          ENDDO
     1602       ENDDO
     1603       !-transform the aerosol forcings, if they have
     1604       ! to be calculated
     1605       IF (ok_ade) THEN
     1606          DO i = 1, kdlon
     1607             topswad_aero(iof+i) = ztopswadaero(i)
     1608             topswad0_aero(iof+i) = ztopswad0aero(i)
     1609             solswad_aero(iof+i) = zsolswadaero(i)
     1610             solswad0_aero(iof+i) = zsolswad0aero(i)
     1611             topsw_aero(iof+i,:) = ztopsw_aero(i,:)
     1612             topsw0_aero(iof+i,:) = ztopsw0_aero(i,:)
     1613             solsw_aero(iof+i,:) = zsolsw_aero(i,:)
     1614             solsw0_aero(iof+i,:) = zsolsw0_aero(i,:)
     1615             topswcf_aero(iof+i,:) = ztopswcf_aero(i,:)
     1616             solswcf_aero(iof+i,:) = zsolswcf_aero(i,:)   
     1617             !-LW
     1618             toplwad_aero(iof+i) = ztoplwadaero(i)
     1619             toplwad0_aero(iof+i) = ztoplwad0aero(i)
     1620             sollwad_aero(iof+i) = zsollwadaero(i)
     1621             sollwad0_aero(iof+i) = zsollwad0aero(i)   
     1622          ENDDO
     1623       ELSE
     1624          DO i = 1, kdlon
     1625             topswad_aero(iof+i) = 0.0
     1626             solswad_aero(iof+i) = 0.0
     1627             topswad0_aero(iof+i) = 0.0
     1628             solswad0_aero(iof+i) = 0.0
     1629             topsw_aero(iof+i,:) = 0.
     1630             topsw0_aero(iof+i,:) =0.
     1631             solsw_aero(iof+i,:) = 0.
     1632             solsw0_aero(iof+i,:) = 0.
     1633             !-LW
     1634             toplwad_aero(iof+i) = 0.0
     1635             sollwad_aero(iof+i) = 0.0
     1636             toplwad0_aero(iof+i) = 0.0
     1637             sollwad0_aero(iof+i) = 0.0
     1638          ENDDO
    15381639       ENDIF
    1539       ENDDO
    1540 !
    1541       DO i = 1, kdlon
    1542          zsolsw(i)    = ZSWFT(i,1)
    1543          zsolsw0(i)   = ZSWFT0_i(i,1)
    1544          ztopsw(i)    = ZSWFT(i,klev+1)
    1545          ztopsw0(i)   = ZSWFT0_i(i,klev+1)
    1546          zsollw(i)    = ZLWFT(i,1)
    1547          zsollw0(i)   = ZLWFT0_i(i,1)
    1548          ztoplw(i)    = ZLWFT(i,klev+1)*(-1)
    1549          ztoplw0(i)   = ZLWFT0_i(i,klev+1)*(-1)
    1550 !         
    1551          zsollwdown(i)= -1.*ZFLDN(i,1)
    1552       ENDDO
    1553 
    1554       DO k=1,kflev
    1555          DO i=1,kdlon
    1556            zheat(i,k)=(ZSWFT(i,k+1)-ZSWFT(i,k))*RDAY*RG/RCPD/PDP(i,k)
    1557            zheat0(i,k)=(ZSWFT0_i(i,k+1)-ZSWFT0_i(i,k))*RDAY*RG/RCPD/PDP(i,k)
    1558            zcool(i,k)=(ZLWFT(i,k)-ZLWFT(i,k+1))*RDAY*RG/RCPD/PDP(i,k)
    1559            zcool0(i,k)=(ZLWFT0_i(i,k)-ZLWFT0_i(i,k+1))*RDAY*RG/RCPD/PDP(i,k)
    1560            IF (ok_volcan) THEN
    1561               zheat_volc(i,k)=(ZSWADAERO(i,k+1)-ZSWADAERO(i,k))*RG/RCPD/PDP(i,k) !NL
    1562               zcool_volc(i,k)=(ZLWADAERO(i,k)-ZLWADAERO(i,k+1))*RG/RCPD/PDP(i,k) !NL
    1563            ENDIF
    1564          ENDDO
    1565       ENDDO
    1566 #endif 
    1567   print*,'Fin traitement ECRAD'
    1568 ! Fin ECRAD
    1569   ENDIF        ! iflag_rrtm
    1570 ! ecrad
    1571 !======================================================================
    1572 
    1573     DO i = 1, kdlon
    1574       topsw(iof+i) = ztopsw(i)
    1575       toplw(iof+i) = ztoplw(i)
    1576       solsw(iof+i) = zsolsw(i)
    1577       solswfdiff(iof+i) = zsolswfdiff(i)
    1578       sollw(iof+i) = zsollw(i)
    1579       sollwdown(iof+i) = zsollwdown(i)
    1580       DO k = 1, kflev+1
    1581         lwdn0 ( iof+i,k)   = ZFLDN0 ( i,k)
    1582         lwdn  ( iof+i,k)   = ZFLDN  ( i,k)
    1583         lwup0 ( iof+i,k)   = ZFLUP0 ( i,k)
    1584         lwup  ( iof+i,k)   = ZFLUP  ( i,k)
    1585       ENDDO
    1586       topsw0(iof+i) = ztopsw0(i)
    1587       toplw0(iof+i) = ztoplw0(i)
    1588       solsw0(iof+i) = zsolsw0(i)
    1589       sollw0(iof+i) = zsollw0(i)
    1590       albpla(iof+i) = zalbpla(i)
    1591 
    1592       DO k = 1, kflev+1
    1593         swdnc0( iof+i,k)   = ZFSDNC0( i,k)
    1594         swdn0 ( iof+i,k)   = ZFSDN0 ( i,k)
    1595         swdn  ( iof+i,k)   = ZFSDN  ( i,k)
    1596         swupc0( iof+i,k)   = ZFSUPC0( i,k)
    1597         swup0 ( iof+i,k)   = ZFSUP0 ( i,k)
    1598         swup  ( iof+i,k)   = ZFSUP  ( i,k)
    1599         lwdnc0( iof+i,k)   = ZFLDNC0( i,k)
    1600         lwupc0( iof+i,k)   = ZFLUPC0( i,k)
    1601       ENDDO
    1602     ENDDO
    1603     !-transform the aerosol forcings, if they have
    1604     ! to be calculated
    1605     IF (ok_ade) THEN
    1606         DO i = 1, kdlon
    1607           topswad_aero(iof+i) = ztopswadaero(i)
    1608           topswad0_aero(iof+i) = ztopswad0aero(i)
    1609           solswad_aero(iof+i) = zsolswadaero(i)
    1610           solswad0_aero(iof+i) = zsolswad0aero(i)
    1611           topsw_aero(iof+i,:) = ztopsw_aero(i,:)
    1612           topsw0_aero(iof+i,:) = ztopsw0_aero(i,:)
    1613           solsw_aero(iof+i,:) = zsolsw_aero(i,:)
    1614           solsw0_aero(iof+i,:) = zsolsw0_aero(i,:)
    1615           topswcf_aero(iof+i,:) = ztopswcf_aero(i,:)
    1616           solswcf_aero(iof+i,:) = zsolswcf_aero(i,:)   
    1617           !-LW
    1618           toplwad_aero(iof+i) = ztoplwadaero(i)
    1619           toplwad0_aero(iof+i) = ztoplwad0aero(i)
    1620           sollwad_aero(iof+i) = zsollwadaero(i)
    1621           sollwad0_aero(iof+i) = zsollwad0aero(i)   
    1622         ENDDO
    1623     ELSE
    1624         DO i = 1, kdlon
    1625           topswad_aero(iof+i) = 0.0
    1626           solswad_aero(iof+i) = 0.0
    1627           topswad0_aero(iof+i) = 0.0
    1628           solswad0_aero(iof+i) = 0.0
    1629           topsw_aero(iof+i,:) = 0.
    1630           topsw0_aero(iof+i,:) =0.
    1631           solsw_aero(iof+i,:) = 0.
    1632           solsw0_aero(iof+i,:) = 0.
    1633           !-LW
    1634           toplwad_aero(iof+i) = 0.0
    1635           sollwad_aero(iof+i) = 0.0
    1636           toplwad0_aero(iof+i) = 0.0
    1637           sollwad0_aero(iof+i) = 0.0
    1638         ENDDO
     1640       IF (ok_aie) THEN
     1641          DO i = 1, kdlon
     1642             topswai_aero(iof+i) = ztopswaiaero(i)
     1643             solswai_aero(iof+i) = zsolswaiaero(i)
     1644             !-LW
     1645             toplwai_aero(iof+i) = ztoplwaiaero(i)
     1646             sollwai_aero(iof+i) = zsollwaiaero(i)
     1647          ENDDO
     1648       ELSE
     1649          DO i = 1, kdlon
     1650             topswai_aero(iof+i) = 0.0
     1651             solswai_aero(iof+i) = 0.0
     1652             !-LW
     1653             toplwai_aero(iof+i) = 0.0
     1654             sollwai_aero(iof+i) = 0.0
     1655          ENDDO
     1656       ENDIF
     1657       DO k = 1, kflev
     1658          DO i = 1, kdlon
     1659             !        scale factor to take into account the difference between
     1660             !        dry air and watter vapour scpecifi! heat capacity
     1661             zznormcp=1.0+RVTMP2*PWV(i,k)
     1662             heat(iof+i,k) = zheat(i,k)/zznormcp
     1663             cool(iof+i,k) = zcool(i,k)/zznormcp
     1664             heat0(iof+i,k) = zheat0(i,k)/zznormcp
     1665             cool0(iof+i,k) = zcool0(i,k)/zznormcp
     1666             IF(ok_volcan) THEN !NL
     1667                heat_volc(iof+i,k) = zheat_volc(i,k)/zznormcp
     1668                cool_volc(iof+i,k) = zcool_volc(i,k)/zznormcp
     1669             ENDIF
     1670          ENDDO
     1671       ENDDO
     1672
     1673    ENDDO ! j = 1, nb_gr
     1674
     1675    IF (lldebug) THEN
     1676       if (0.eq.1) then
     1677          ! Verifs dans le cas 1D
     1678          print*,'================== Sortie de radlw ================='
     1679          print*,'******** LW LW LW *******************'
     1680          print*,'ZLWFT =',ZLWFT
     1681          print*,'ZLWFT0_i =',ZLWFT0_i
     1682          print*,'ZFLUP0 =',ZFLUP0
     1683          print*,'ZFLDN0 =',ZFLDN0
     1684          print*,'ZFLDNC0 =',ZFLDNC0
     1685          print*,'ZFLUPC0 =',ZFLUPC0
     1686
     1687          print*,'******** SW SW SW *******************'
     1688          print*,'ZSWFT =',ZSWFT
     1689          print*,'ZSWFT0_i =',ZSWFT0_i
     1690          print*,'ZFSDN =',ZFSDN
     1691          print*,'ZFSDN0 =',ZFSDN0
     1692          print*,'ZFSDNC0 =',ZFSDNC0
     1693          print*,'ZFSUP =',ZFSUP
     1694          print*,'ZFSUP0 =',ZFSUP0
     1695          print*,'ZFSUPC0 =',ZFSUPC0
     1696
     1697          print*,'******** LMDZ  *******************'
     1698          print*,'cool = ', cool
     1699          print*,'heat = ', heat
     1700          print*,'topsw = ', topsw
     1701          print*,'toplw = ', toplw
     1702          print*,'sollw = ', sollw
     1703          print*,'solsw = ', solsw
     1704          print*,'lwdn = ', lwdn
     1705          print*,'lwup = ', lwup
     1706          print*,'swdn = ', swdn
     1707          print*,'swup =', swup
     1708       endif
    16391709    ENDIF
    1640     IF (ok_aie) THEN
    1641         DO i = 1, kdlon
    1642           topswai_aero(iof+i) = ztopswaiaero(i)
    1643           solswai_aero(iof+i) = zsolswaiaero(i)
    1644           !-LW
    1645           toplwai_aero(iof+i) = ztoplwaiaero(i)
    1646           sollwai_aero(iof+i) = zsollwaiaero(i)
    1647         ENDDO
    1648     ELSE
    1649         DO i = 1, kdlon
    1650           topswai_aero(iof+i) = 0.0
    1651           solswai_aero(iof+i) = 0.0
    1652           !-LW
    1653           toplwai_aero(iof+i) = 0.0
    1654           sollwai_aero(iof+i) = 0.0
    1655         ENDDO
    1656     ENDIF
    1657     DO k = 1, kflev
    1658       DO i = 1, kdlon
    1659         !        scale factor to take into account the difference between
    1660         !        dry air and watter vapour scpecifi! heat capacity
    1661         zznormcp=1.0+RVTMP2*PWV(i,k)
    1662         heat(iof+i,k) = zheat(i,k)/zznormcp
    1663         cool(iof+i,k) = zcool(i,k)/zznormcp
    1664         heat0(iof+i,k) = zheat0(i,k)/zznormcp
    1665         cool0(iof+i,k) = zcool0(i,k)/zznormcp
    1666         IF(ok_volcan) THEN !NL
    1667            heat_volc(iof+i,k) = zheat_volc(i,k)/zznormcp
    1668            cool_volc(iof+i,k) = zcool_volc(i,k)/zznormcp
    1669         ENDIF
    1670       ENDDO
    1671     ENDDO
    1672 
    1673  ENDDO ! j = 1, nb_gr
    1674 
    1675 IF (lldebug) THEN
    1676  if (0.eq.1) then
    1677 ! Verifs dans le cas 1D
    1678  print*,'================== Sortie de radlw ================='
    1679  print*,'******** LW LW LW *******************'
    1680  print*,'ZLWFT =',ZLWFT
    1681  print*,'ZLWFT0_i =',ZLWFT0_i
    1682  print*,'ZFLUP0 =',ZFLUP0
    1683  print*,'ZFLDN0 =',ZFLDN0
    1684  print*,'ZFLDNC0 =',ZFLDNC0
    1685  print*,'ZFLUPC0 =',ZFLUPC0
    1686 
    1687  print*,'******** SW SW SW *******************'
    1688  print*,'ZSWFT =',ZSWFT
    1689  print*,'ZSWFT0_i =',ZSWFT0_i
    1690  print*,'ZFSDN =',ZFSDN
    1691  print*,'ZFSDN0 =',ZFSDN0
    1692  print*,'ZFSDNC0 =',ZFSDNC0
    1693  print*,'ZFSUP =',ZFSUP
    1694  print*,'ZFSUP0 =',ZFSUP0
    1695  print*,'ZFSUPC0 =',ZFSUPC0
    1696 
    1697  print*,'******** LMDZ  *******************'
    1698  print*,'cool = ', cool
    1699  print*,'heat = ', heat
    1700  print*,'topsw = ', topsw
    1701  print*,'toplw = ', toplw
    1702  print*,'sollw = ', sollw
    1703  print*,'solsw = ', solsw
    1704  print*,'lwdn = ', lwdn
    1705  print*,'lwup = ', lwup
    1706  print*,'swdn = ', swdn
    1707  print*,'swup =', swup
    1708  endif
    1709 ENDIF
    1710 
    1711 END SUBROUTINE radlwsw
     1710
     1711  END SUBROUTINE radlwsw
    17121712
    17131713end module radlwsw_m
Note: See TracChangeset for help on using the changeset viewer.