Context Navigation

radlwsw_m.F90 @ 4031

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

Corrections in the LMDZ-ECRAN interface:

gas concentrations in VMR instead of MMR
use of a constant value for each gas instead of a table of values
specify in the 14 SW bands of ECRAD those corresponding to the 6 SW bands for diffuse and direct albedo in LMDZ
passage of all processing of input fields for Ecrad (conversions, calculations, ...) from radiation.F90 to radlw.F90

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:keywords set to Author Date Id Revi

File size: 62.7 KB

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

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

Download in other formats:

Original Format