source: LMDZ6/trunk/libf/phylmd/radlwsw_m.F90 @ 4046

Last change on this file since 4046 was 4045, checked in by idelkadi, 3 years ago

Correction in radlwsw_m.F90:
Addition of machine threshold to avoid a division by 0

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