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

Last change on this file since 4853 was 4853, checked in by idelkadi, 2 months ago

Ecrad update in LMDZ / Implementation of Ecrad double call in LMDZ

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