Context Navigation

radlwsw_m.F90 @ 4754

Last change on this file since 4754 was 4677, checked in by idelkadi, 15 months ago

Implementation in the LMDZ code of the double call of the ECRAD radiative transfer code to estimate the 3D radiative effect of clouds.

This double call of Ecrad is controlled by the ok_3Deffect logic key.
If this key is enabled, 2 files of parameter configuration "namelists" for ECRAD are required at runtime: namelist_ecrad and namelist_ecrad_s2.
If this key is deactivated, the configuration and initialization part (reading namelist and netcdf files) is performed only once during simulation (1st call to ECRAD). Otherwise, configuration and initialization are performed each time Ecrad is called.

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.6 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format