Context Navigation

radlwsw_m.F90 @ 4283

Last change on this file since 4283 was 4170, checked in by dcugnet, 2 years ago
The variable "types_trac" is the equivalent of "type_trac" in case multiple sections must be read and used in "tracer.def" file. Tests on the "type_trac" were replaced with tests on the vector "types_trac". Most of the time, there are two components: 'lmdz' and a second one. The later has priority on 'lmdz' and must be used for the tests. For more components, care must be taken to execute specific parts of the code on the right tracers ; the tracers(:)%component has been created in that respect.
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:keywords set to `Author Date Id Revi`
File size: 63.5 KB

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

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

Download in other formats:

Original Format