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

Last change on this file since 4163 was 4116, checked in by idelkadi, 3 years ago

Implementation of Ecrad in LMDZ (continued):

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