Context Navigation

radlwsw_m.F90 @ 3985

Last change on this file since 3985 was 3406, checked in by jghattas, 6 years ago
Added all modifications in the model code that were used for the simulations with DYANMICO during the Grand Challeng 2018. Modifications done by Y. Meurdesoif, L. Fairhead and A.K. Traore
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: 44.3 KB

Line
1	!
2	! $Id: radlwsw_m.F90 3406 2018-10-23 09:16:41Z dcugnet $
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, flag_aerosol,&
19	flag_aerosol_strat,&
20	tau_aero, piz_aero, cg_aero,&
21	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm,& ! rajoute par OB pour RRTM
22	tau_aero_lw_rrtm, & ! rajoute par C. Kleinschmitt pour RRTM
23	cldtaupi, new_aod, &
24	qsat, flwc, fiwc, &
25	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
26	heat,heat0,cool,cool0,albpla,&
27	topsw,toplw,solsw,sollw,&
28	sollwdown,&
29	topsw0,toplw0,solsw0,sollw0,&
30	lwdn0, lwdn, lwup0, lwup,&
31	swdn0, swdn, swup0, swup,&
32	topswad_aero, solswad_aero,&
33	topswai_aero, solswai_aero, &
34	topswad0_aero, solswad0_aero,&
35	topsw_aero, topsw0_aero,&
36	solsw_aero, solsw0_aero, &
37	topswcf_aero, solswcf_aero,&
38	!-C. Kleinschmitt for LW diagnostics
39	toplwad_aero, sollwad_aero,&
40	toplwai_aero, sollwai_aero, &
41	toplwad0_aero, sollwad0_aero,&
42	!-end
43	ZLWFT0_i, ZFLDN0, ZFLUP0,&
44	ZSWFT0_i, ZFSDN0, ZFSUP0)
45
46
47
48	USE DIMPHY
49	USE assert_m, ONLY : assert
50	USE infotrac_phy, ONLY : type_trac
51	USE write_field_phy
52	#ifdef REPROBUS
53	USE CHEM_REP, ONLY : solaireTIME, ok_SUNTIME, ndimozon
54	#endif
55	#ifdef CPP_RRTM
56	! modules necessaires au rayonnement
57	! -----------------------------------------
58	! USE YOMCST , ONLY : RG ,RD ,RTT ,RPI
59	! USE YOERAD , ONLY : NSW ,LRRTM ,LINHOM , LCCNL,LCCNO,
60	! USE YOERAD , ONLY : NSW ,LRRTM ,LCCNL ,LCCNO ,&
61	! NSW mis dans .def MPL 20140211
62	! NLW ajoute par OB
63	USE YOERAD , ONLY : NLW, LRRTM ,LCCNL ,LCCNO ,&
64	NRADIP , NRADLP , NICEOPT, NLIQOPT ,RCCNLND , RCCNSEA
65	USE YOELW , ONLY : NSIL ,NTRA ,NUA ,TSTAND ,XP
66	USE YOESW , ONLY : RYFWCA ,RYFWCB ,RYFWCC ,RYFWCD,&
67	RYFWCE ,RYFWCF ,REBCUA ,REBCUB ,REBCUC,&
68	REBCUD ,REBCUE ,REBCUF ,REBCUI ,REBCUJ,&
69	REBCUG ,REBCUH ,RHSAVI ,RFULIO ,RFLAA0,&
70	RFLAA1 ,RFLBB0 ,RFLBB1 ,RFLBB2 ,RFLBB3,&
71	RFLCC0 ,RFLCC1 ,RFLCC2 ,RFLCC3 ,RFLDD0,&
72	RFLDD1 ,RFLDD2 ,RFLDD3 ,RFUETA ,RASWCA,&
73	RASWCB ,RASWCC ,RASWCD ,RASWCE ,RASWCF
74	! & RASWCB ,RASWCC ,RASWCD ,RASWCE ,RASWCF, RLINLI
75	USE YOERDU , ONLY : NUAER ,NTRAER ,REPLOG ,REPSC ,REPSCW ,DIFF
76	! USE YOETHF , ONLY : RTICE
77	USE YOERRTWN , ONLY : DELWAVE ,TOTPLNK
78	USE YOMPHY3 , ONLY : RII0
79	#endif
80	USE aero_mod
81
82	!======================================================================
83	! Auteur(s): Z.X. Li (LMD/CNRS) date: 19960719
84	! Objet: interface entre le modele et les rayonnements
85	! Arguments:
86	! dist-----input-R- distance astronomique terre-soleil
87	! rmu0-----input-R- cosinus de l'angle zenithal
88	! fract----input-R- duree d'ensoleillement normalisee
89	! co2_ppm--input-R- concentration du gaz carbonique (en ppm)
90	! paprs----input-R- pression a inter-couche (Pa)
91	! pplay----input-R- pression au milieu de couche (Pa)
92	! tsol-----input-R- temperature du sol (en K)
93	! alb1-----input-R- albedo du sol(entre 0 et 1) dans l'interval visible
94	! alb2-----input-R- albedo du sol(entre 0 et 1) dans l'interval proche infra-rouge
95	! t--------input-R- temperature (K)
96	! q--------input-R- vapeur d'eau (en kg/kg)
97	! cldfra---input-R- fraction nuageuse (entre 0 et 1)
98	! cldtaupd---input-R- epaisseur optique des nuages dans le visible (present-day value)
99	! cldemi---input-R- emissivite des nuages dans l'IR (entre 0 et 1)
100	! ok_ade---input-L- apply the Aerosol Direct Effect or not?
101	! ok_aie---input-L- apply the Aerosol Indirect Effect or not?
102	! flag_aerosol-input-I- aerosol flag from 0 to 6
103	! flag_aerosol_strat-input-I- use stratospheric aerosols flag (0, 1, 2)
104	! tau_ae, piz_ae, cg_ae-input-R- aerosol optical properties (calculated in aeropt.F)
105	! cldtaupi-input-R- epaisseur optique des nuages dans le visible
106	! calculated for pre-industrial (pi) aerosol concentrations, i.e. with smaller
107	! droplet concentration, thus larger droplets, thus generally cdltaupi cldtaupd
108	! it is needed for the diagnostics of the aerosol indirect radiative forcing
109	!
110	! heat-----output-R- echauffement atmospherique (visible) (K/jour)
111	! cool-----output-R- refroidissement dans l'IR (K/jour)
112	! albpla---output-R- albedo planetaire (entre 0 et 1)
113	! topsw----output-R- flux solaire net au sommet de l'atm.
114	! toplw----output-R- ray. IR montant au sommet de l'atmosphere
115	! solsw----output-R- flux solaire net a la surface
116	! sollw----output-R- ray. IR montant a la surface
117	! solswad---output-R- ray. solaire net absorbe a la surface (aerosol dir)
118	! topswad---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol dir)
119	! solswai---output-R- ray. solaire net absorbe a la surface (aerosol ind)
120	! topswai---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol ind)
121	!
122	! ATTENTION: swai and swad have to be interpreted in the following manner:
123	! ---------
124	! ok_ade=F & ok_aie=F -both are zero
125	! ok_ade=T & ok_aie=F -aerosol direct forcing is F_{AD} = topsw-topswad
126	! indirect is zero
127	! ok_ade=F & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai
128	! direct is zero
129	! ok_ade=T & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai
130	! aerosol direct forcing is F_{AD} = topswai-topswad
131	!
132	! --------- RRTM: output RECMWFL
133	! ZEMTD (KPROMA,KLEV+1) ; TOTAL DOWNWARD LONGWAVE EMISSIVITY
134	! ZEMTU (KPROMA,KLEV+1) ; TOTAL UPWARD LONGWAVE EMISSIVITY
135	! ZTRSO (KPROMA,KLEV+1) ; TOTAL SHORTWAVE TRANSMISSIVITY
136	! ZTH (KPROMA,KLEV+1) ; HALF LEVEL TEMPERATURE
137	! ZCTRSO(KPROMA,2) ; CLEAR-SKY SHORTWAVE TRANSMISSIVITY
138	! ZCEMTR(KPROMA,2) ; CLEAR-SKY NET LONGWAVE EMISSIVITY
139	! ZTRSOD(KPROMA) ; TOTAL-SKY SURFACE SW TRANSMISSITY
140	! ZLWFC (KPROMA,2) ; CLEAR-SKY LONGWAVE FLUXES
141	! ZLWFT (KPROMA,KLEV+1) ; TOTAL-SKY LONGWAVE FLUXES
142	! ZLWFT0(KPROMA,KLEV+1) ; CLEAR-SKY LONGWAVE FLUXES ! added by MPL 090109
143	! ZSWFC (KPROMA,2) ; CLEAR-SKY SHORTWAVE FLUXES
144	! ZSWFT (KPROMA,KLEV+1) ; TOTAL-SKY SHORTWAVE FLUXES
145	! ZSWFT0(KPROMA,KLEV+1) ; CLEAR-SKY SHORTWAVE FLUXES ! added by MPL 090109
146	! ZFLUX (KLON,2,KLEV+1) ; TOTAL LW FLUXES 1=up, 2=DWN ! added by MPL 080411
147	! ZFLUC (KLON,2,KLEV+1) ; CLEAR SKY LW FLUXES ! added by MPL 080411
148	! ZFSDWN(klon,KLEV+1) ; TOTAL SW DWN FLUXES ! added by MPL 080411
149	! ZFCDWN(klon,KLEV+1) ; CLEAR SKY SW DWN FLUXES ! added by MPL 080411
150	! ZFSUP (klon,KLEV+1) ; TOTAL SW UP FLUXES ! added by MPL 080411
151	! ZFCUP (klon,KLEV+1) ; CLEAR SKY SW UP FLUXES ! added by MPL 080411
152
153	!======================================================================
154
155	! ====================================================================
156	! Adapte au modele de chimie INCA par Celine Deandreis & Anne Cozic -- 2009
157	! 1 = ZERO
158	! 2 = AER total
159	! 3 = NAT
160	! 4 = BC
161	! 5 = SO4
162	! 6 = POM
163	! 7 = DUST
164	! 8 = SS
165	! 9 = NO3
166	!
167	! ====================================================================
168	include "YOETHF.h"
169	include "YOMCST.h"
170	include "clesphys.h"
171
172	! Input arguments
173	REAL, INTENT(in) :: dist
174	REAL, INTENT(in) :: rmu0(KLON), fract(KLON)
175	REAL, INTENT(in) :: paprs(KLON,KLEV+1), pplay(KLON,KLEV)
176	!albedo SB >>>
177	! REAL, INTENT(in) :: alb1(KLON), alb2(KLON), tsol(KLON)
178	REAL, INTENT(in) :: tsol(KLON)
179	REAL, INTENT(in) :: alb_dir(KLON,NSW),alb_dif(KLON,NSW)
180	real, intent(in) :: SFRWL(6)
181	!albedo SB <<<
182	REAL, INTENT(in) :: t(KLON,KLEV), q(KLON,KLEV)
183
184	REAL, INTENT(in):: wo(:, :, :) ! dimension(KLON,KLEV, 1 or 2)
185	! column-density of ozone in a layer, in kilo-Dobsons
186	! "wo(:, :, 1)" is for the average day-night field,
187	! "wo(:, :, 2)" is for daylight time.
188
189	LOGICAL, INTENT(in) :: ok_ade, ok_aie ! switches whether to use aerosol direct (indirect) effects or not
190	LOGICAL :: lldebug
191	INTEGER, INTENT(in) :: flag_aerosol ! takes value 0 (no aerosol) or 1 to 6 (aerosols)
192	INTEGER, INTENT(in) :: flag_aerosol_strat ! use stratospheric aerosols
193	REAL, INTENT(in) :: cldfra(KLON,KLEV), cldemi(KLON,KLEV), cldtaupd(KLON,KLEV)
194	REAL, INTENT(in) :: tau_aero(KLON,KLEV,naero_grp,2) ! aerosol optical properties (see aeropt.F)
195	REAL, INTENT(in) :: piz_aero(KLON,KLEV,naero_grp,2) ! aerosol optical properties (see aeropt.F)
196	REAL, INTENT(in) :: cg_aero(KLON,KLEV,naero_grp,2) ! aerosol optical properties (see aeropt.F)
197	!--OB
198	REAL, INTENT(in) :: tau_aero_sw_rrtm(KLON,KLEV,2,NSW) ! aerosol optical properties RRTM
199	REAL, INTENT(in) :: piz_aero_sw_rrtm(KLON,KLEV,2,NSW) ! aerosol optical properties RRTM
200	REAL, INTENT(in) :: cg_aero_sw_rrtm(KLON,KLEV,2,NSW) ! aerosol optical properties RRTM
201	!--OB fin
202
203	!--C. Kleinschmitt
204	#ifdef CPP_RRTM
205	REAL, INTENT(in) :: tau_aero_lw_rrtm(KLON,KLEV,2,NLW) ! LW aerosol optical properties RRTM
206	#else
207	REAL, INTENT(in) :: tau_aero_lw_rrtm(KLON,KLEV,2,nbands_lw_rrtm)
208	#endif
209	!--C. Kleinschmitt end
210
211	REAL, INTENT(in) :: cldtaupi(KLON,KLEV) ! cloud optical thickness for pre-industrial aerosol concentrations
212	LOGICAL, INTENT(in) :: new_aod ! flag pour retrouver les resultats exacts de l'AR4 dans le cas ou l'on ne travaille qu'avec les sulfates
213	REAL, INTENT(in) :: qsat(klon,klev) ! Variable pour iflag_rrtm=1
214	REAL, INTENT(in) :: flwc(klon,klev) ! Variable pour iflag_rrtm=1
215	REAL, INTENT(in) :: fiwc(klon,klev) ! Variable pour iflag_rrtm=1
216	REAL, INTENT(in) :: ref_liq(klon,klev) ! cloud droplet radius present-day from newmicro
217	REAL, INTENT(in) :: ref_ice(klon,klev) ! ice crystal radius present-day from newmicro
218	REAL, INTENT(in) :: ref_liq_pi(klon,klev) ! cloud droplet radius pre-industrial from newmicro
219	REAL, INTENT(in) :: ref_ice_pi(klon,klev) ! ice crystal radius pre-industrial from newmicro
220
221	! Output arguments
222	REAL, INTENT(out) :: heat(KLON,KLEV), cool(KLON,KLEV)
223	REAL, INTENT(out) :: heat0(KLON,KLEV), cool0(KLON,KLEV)
224	REAL, INTENT(out) :: topsw(KLON), toplw(KLON)
225	REAL, INTENT(out) :: solsw(KLON), sollw(KLON), albpla(KLON)
226	REAL, INTENT(out) :: topsw0(KLON), toplw0(KLON), solsw0(KLON), sollw0(KLON)
227	REAL, INTENT(out) :: sollwdown(KLON)
228	REAL, INTENT(out) :: swdn(KLON,kflev+1),swdn0(KLON,kflev+1)
229	REAL, INTENT(out) :: swup(KLON,kflev+1),swup0(KLON,kflev+1)
230	REAL, INTENT(out) :: lwdn(KLON,kflev+1),lwdn0(KLON,kflev+1)
231	REAL, INTENT(out) :: lwup(KLON,kflev+1),lwup0(KLON,kflev+1)
232	REAL, INTENT(out) :: topswad_aero(KLON), solswad_aero(KLON) ! output: aerosol direct forcing at TOA and surface
233	REAL, INTENT(out) :: topswai_aero(KLON), solswai_aero(KLON) ! output: aerosol indirect forcing atTOA and surface
234	REAL, INTENT(out) :: toplwad_aero(KLON), sollwad_aero(KLON) ! output: LW aerosol direct forcing at TOA and surface
235	REAL, INTENT(out) :: toplwai_aero(KLON), sollwai_aero(KLON) ! output: LW aerosol indirect forcing atTOA and surface
236	REAL, DIMENSION(klon), INTENT(out) :: topswad0_aero
237	REAL, DIMENSION(klon), INTENT(out) :: solswad0_aero
238	REAL, DIMENSION(klon), INTENT(out) :: toplwad0_aero
239	REAL, DIMENSION(klon), INTENT(out) :: sollwad0_aero
240	REAL, DIMENSION(kdlon,9), INTENT(out) :: topsw_aero
241	REAL, DIMENSION(kdlon,9), INTENT(out) :: topsw0_aero
242	REAL, DIMENSION(kdlon,9), INTENT(out) :: solsw_aero
243	REAL, DIMENSION(kdlon,9), INTENT(out) :: solsw0_aero
244	REAL, DIMENSION(kdlon,3), INTENT(out) :: topswcf_aero
245	REAL, DIMENSION(kdlon,3), INTENT(out) :: solswcf_aero
246	REAL, DIMENSION(kdlon,kflev+1), INTENT(out) :: ZSWFT0_i
247	REAL, DIMENSION(kdlon,kflev+1), INTENT(out) :: ZLWFT0_i
248
249	! Local variables
250	REAL(KIND=8) ZFSUP(KDLON,KFLEV+1)
251	REAL(KIND=8) ZFSDN(KDLON,KFLEV+1)
252	REAL(KIND=8) ZFSUP0(KDLON,KFLEV+1)
253	REAL(KIND=8) ZFSDN0(KDLON,KFLEV+1)
254	REAL(KIND=8) ZFLUP(KDLON,KFLEV+1)
255	REAL(KIND=8) ZFLDN(KDLON,KFLEV+1)
256	REAL(KIND=8) ZFLUP0(KDLON,KFLEV+1)
257	REAL(KIND=8) ZFLDN0(KDLON,KFLEV+1)
258	REAL(KIND=8) zx_alpha1, zx_alpha2
259	INTEGER k, kk, i, j, iof, nb_gr
260	INTEGER ist,iend,ktdia,kmode
261	REAL(KIND=8) PSCT
262	REAL(KIND=8) PALBD(kdlon,2), PALBP(kdlon,2)
263	! MPL 06.01.09: pour RRTM, creation de PALBD_NEW et PALBP_NEW
264	! avec NSW en deuxieme dimension
265	REAL(KIND=8) PALBD_NEW(kdlon,NSW), PALBP_NEW(kdlon,NSW)
266	REAL(KIND=8) PEMIS(kdlon), PDT0(kdlon), PVIEW(kdlon)
267	REAL(KIND=8) PPSOL(kdlon), PDP(kdlon,KLEV)
268	REAL(KIND=8) PTL(kdlon,kflev+1), PPMB(kdlon,kflev+1)
269	REAL(KIND=8) PTAVE(kdlon,kflev)
270	REAL(KIND=8) PWV(kdlon,kflev), PQS(kdlon,kflev)
271
272	real(kind=8) POZON(kdlon, kflev, size(wo, 3)) ! mass fraction of ozone
273	! "POZON(:, :, 1)" is for the average day-night field,
274	! "POZON(:, :, 2)" is for daylight time.
275	!!!!! Modif MPL 6.01.09 avec RRTM, on passe de 5 a 6
276	REAL(KIND=8) PAER(kdlon,kflev,6)
277	REAL(KIND=8) PCLDLD(kdlon,kflev)
278	REAL(KIND=8) PCLDLU(kdlon,kflev)
279	REAL(KIND=8) PCLDSW(kdlon,kflev)
280	REAL(KIND=8) PTAU(kdlon,2,kflev)
281	REAL(KIND=8) POMEGA(kdlon,2,kflev)
282	REAL(KIND=8) PCG(kdlon,2,kflev)
283	REAL(KIND=8) zfract(kdlon), zrmu0(kdlon), zdist
284	REAL(KIND=8) zheat(kdlon,kflev), zcool(kdlon,kflev)
285	REAL(KIND=8) zheat0(kdlon,kflev), zcool0(kdlon,kflev)
286	REAL(KIND=8) ztopsw(kdlon), ztoplw(kdlon)
287	REAL(KIND=8) zsolsw(kdlon), zsollw(kdlon), zalbpla(kdlon)
288	REAL(KIND=8) zsollwdown(kdlon)
289	REAL(KIND=8) ztopsw0(kdlon), ztoplw0(kdlon)
290	REAL(KIND=8) zsolsw0(kdlon), zsollw0(kdlon)
291	REAL(KIND=8) zznormcp
292	REAL(KIND=8) tauaero(kdlon,kflev,naero_grp,2) ! aer opt properties
293	REAL(KIND=8) pizaero(kdlon,kflev,naero_grp,2)
294	REAL(KIND=8) cgaero(kdlon,kflev,naero_grp,2)
295	REAL(KIND=8) PTAUA(kdlon,2,kflev) ! present-day value of cloud opt thickness (PTAU is pre-industrial value), local use
296	REAL(KIND=8) POMEGAA(kdlon,2,kflev) ! dito for single scatt albedo
297	REAL(KIND=8) ztopswadaero(kdlon), zsolswadaero(kdlon) ! Aerosol direct forcing at TOAand surface
298	REAL(KIND=8) ztopswad0aero(kdlon), zsolswad0aero(kdlon) ! Aerosol direct forcing at TOAand surface
299	REAL(KIND=8) ztopswaiaero(kdlon), zsolswaiaero(kdlon) ! dito, indirect
300	!-LW by CK
301	REAL(KIND=8) ztoplwadaero(kdlon), zsollwadaero(kdlon) ! LW Aerosol direct forcing at TOAand surface
302	REAL(KIND=8) ztoplwad0aero(kdlon), zsollwad0aero(kdlon) ! LW Aerosol direct forcing at TOAand surface
303	REAL(KIND=8) ztoplwaiaero(kdlon), zsollwaiaero(kdlon) ! dito, indirect
304	!-end
305	REAL(KIND=8) ztopsw_aero(kdlon,9), ztopsw0_aero(kdlon,9)
306	REAL(KIND=8) zsolsw_aero(kdlon,9), zsolsw0_aero(kdlon,9)
307	REAL(KIND=8) ztopswcf_aero(kdlon,3), zsolswcf_aero(kdlon,3)
308	! real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2 deje declare dans physiq.F MPL 20130618
309	!MPL input supplementaires pour RECMWFL
310	! flwc, fiwc = Liquid Water Content & Ice Water Content (kg/kg)
311	REAL(KIND=8) GEMU(klon)
312	!MPL input RECMWFL:
313	! Tableaux aux niveaux inverses pour respecter convention Arpege
314	REAL(KIND=8) ref_liq_i(klon,klev) ! cloud droplet radius present-day from newmicro (inverted)
315	REAL(KIND=8) ref_ice_i(klon,klev) ! ice crystal radius present-day from newmicro (inverted)
316	!--OB
317	REAL(KIND=8) ref_liq_pi_i(klon,klev) ! cloud droplet radius pre-industrial from newmicro (inverted)
318	REAL(KIND=8) ref_ice_pi_i(klon,klev) ! ice crystal radius pre-industrial from newmicro (inverted)
319	!--end OB
320	REAL(KIND=8) paprs_i(klon,klev+1)
321	REAL(KIND=8) pplay_i(klon,klev)
322	REAL(KIND=8) cldfra_i(klon,klev)
323	REAL(KIND=8) POZON_i(kdlon,kflev, size(wo, 3)) ! mass fraction of ozone
324	! "POZON(:, :, 1)" is for the average day-night field,
325	! "POZON(:, :, 2)" is for daylight time.
326	!!!!! Modif MPL 6.01.09 avec RRTM, on passe de 5 a 6
327	REAL(KIND=8) PAER_i(kdlon,kflev,6)
328	REAL(KIND=8) PDP_i(klon,klev)
329	REAL(KIND=8) t_i(klon,klev),q_i(klon,klev),qsat_i(klon,klev)
330	REAL(KIND=8) flwc_i(klon,klev),fiwc_i(klon,klev)
331	!MPL output RECMWFL:
332	REAL(KIND=8) ZEMTD (klon,klev+1),ZEMTD_i (klon,klev+1)
333	REAL(KIND=8) ZEMTU (klon,klev+1),ZEMTU_i (klon,klev+1)
334	REAL(KIND=8) ZTRSO (klon,klev+1),ZTRSO_i (klon,klev+1)
335	REAL(KIND=8) ZTH (klon,klev+1),ZTH_i (klon,klev+1)
336	REAL(KIND=8) ZCTRSO(klon,2)
337	REAL(KIND=8) ZCEMTR(klon,2)
338	REAL(KIND=8) ZTRSOD(klon)
339	REAL(KIND=8) ZLWFC (klon,2)
340	REAL(KIND=8) ZLWFT (klon,klev+1),ZLWFT_i (klon,klev+1)
341	REAL(KIND=8) ZSWFC (klon,2)
342	REAL(KIND=8) ZSWFT (klon,klev+1),ZSWFT_i (klon,klev+1)
343	REAL(KIND=8) ZFLUCDWN_i(klon,klev+1),ZFLUCUP_i(klon,klev+1)
344	REAL(KIND=8) PPIZA_TOT(klon,klev,NSW)
345	REAL(KIND=8) PCGA_TOT(klon,klev,NSW)
346	REAL(KIND=8) PTAU_TOT(klon,klev,NSW)
347	REAL(KIND=8) PPIZA_NAT(klon,klev,NSW)
348	REAL(KIND=8) PCGA_NAT(klon,klev,NSW)
349	REAL(KIND=8) PTAU_NAT(klon,klev,NSW)
350	#ifdef CPP_RRTM
351	REAL(KIND=8) PTAU_LW_TOT(klon,klev,NLW)
352	REAL(KIND=8) PTAU_LW_NAT(klon,klev,NLW)
353	#endif
354	REAL(KIND=8) PSFSWDIR(klon,NSW)
355	REAL(KIND=8) PSFSWDIF(klon,NSW)
356	REAL(KIND=8) PFSDNN(klon)
357	REAL(KIND=8) PFSDNV(klon)
358	!MPL On ne redefinit pas les tableaux ZFLUX,ZFLUC,
359	!MPL ZFSDWN,ZFCDWN,ZFSUP,ZFCUP car ils existent deja
360	!MPL sous les noms de ZFLDN,ZFLDN0,ZFLUP,ZFLUP0,
361	!MPL ZFSDN,ZFSDN0,ZFSUP,ZFSUP0
362	REAL(KIND=8) ZFLUX_i (klon,2,klev+1)
363	REAL(KIND=8) ZFLUC_i (klon,2,klev+1)
364	REAL(KIND=8) ZFSDWN_i (klon,klev+1)
365	REAL(KIND=8) ZFCDWN_i (klon,klev+1)
366	REAL(KIND=8) ZFSUP_i (klon,klev+1)
367	REAL(KIND=8) ZFCUP_i (klon,klev+1)
368	! 3 lignes suivantes a activer pour CCMVAL (MPL 20100412)
369	! REAL(KIND=8) RSUN(3,2)
370	! REAL(KIND=8) SUN(3)
371	! REAL(KIND=8) SUN_FRACT(2)
372	real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
373	CHARACTER (LEN=80) :: abort_message
374	CHARACTER (LEN=80) :: modname='radlwsw_m'
375
376	call assert(size(wo, 1) == klon, size(wo, 2) == klev, "radlwsw wo")
377	! initialisation
378	ist=1
379	iend=klon
380	ktdia=1
381	kmode=ist
382	tauaero(:,:,:,:)=0.
383	pizaero(:,:,:,:)=0.
384	cgaero(:,:,:,:)=0.
385	lldebug=.FALSE.
386
387	ztopsw_aero(:,:) = 0. !ym missing init : warning : not initialized in SW_AEROAR4
388	ztopsw0_aero(:,:) = 0. !ym missing init : warning : not initialized in SW_AEROAR4
389	zsolsw_aero(:,:) = 0. !ym missing init : warning : not initialized in SW_AEROAR4
390	zsolsw0_aero(:,:) = 0. !ym missing init : warning : not initialized in SW_AEROAR4
391
392	!
393	!-------------------------------------------
394	nb_gr = KLON / kdlon
395	IF (nb_gr*kdlon .NE. KLON) THEN
396	PRINT*, "kdlon mauvais:", KLON, kdlon, nb_gr
397	call abort_physic("radlwsw", "", 1)
398	ENDIF
399	IF (kflev .NE. KLEV) THEN
400	PRINT*, "kflev differe de KLEV, kflev, KLEV"
401	call abort_physic("radlwsw", "", 1)
402	ENDIF
403	!-------------------------------------------
404	DO k = 1, KLEV
405	DO i = 1, KLON
406	heat(i,k)=0.
407	cool(i,k)=0.
408	heat0(i,k)=0.
409	cool0(i,k)=0.
410	ENDDO
411	ENDDO
412	!
413	zdist = dist
414	!
415	PSCT = solaire/zdist/zdist
416
417	IF (type_trac == 'repr') THEN
418	#ifdef REPROBUS
419	if(ok_SUNTIME) PSCT = solaireTIME/zdist/zdist
420	print,'Constante solaire: ',PSCTzdist*zdist
421	#endif
422	END IF
423
424	DO j = 1, nb_gr
425	iof = kdlon*(j-1)
426	DO i = 1, kdlon
427	zfract(i) = fract(iof+i)
428	! zfract(i) = 1. !!!!!! essai MPL 19052010
429	zrmu0(i) = rmu0(iof+i)
430
431
432	!albedo SB >>>
433	!
434	IF (iflag_rrtm==0) THEN
435	!
436	PALBD(i,1)=alb_dif(iof+i,1)
437	PALBD(i,2)=alb_dif(iof+i,2)
438	PALBP(i,1)=alb_dir(iof+i,1)
439	PALBP(i,2)=alb_dir(iof+i,2)
440	!
441	ELSEIF (iflag_rrtm==1) THEn
442	!
443	DO kk=1,NSW
444	PALBD_NEW(i,kk)=alb_dif(iof+i,kk)
445	PALBP_NEW(i,kk)=alb_dir(iof+i,kk)
446	ENDDO
447	!
448	ENDIF
449	!albedo SB <<<
450
451
452	PEMIS(i) = 1.0 !!!!! A REVOIR (MPL)
453	PVIEW(i) = 1.66
454	PPSOL(i) = paprs(iof+i,1)
455	zx_alpha1 = (paprs(iof+i,1)-pplay(iof+i,2))/(pplay(iof+i,1)-pplay(iof+i,2))
456	zx_alpha2 = 1.0 - zx_alpha1
457	PTL(i,1) = t(iof+i,1) * zx_alpha1 + t(iof+i,2) * zx_alpha2
458	PTL(i,KLEV+1) = t(iof+i,KLEV)
459	PDT0(i) = tsol(iof+i) - PTL(i,1)
460	ENDDO
461	DO k = 2, kflev
462	DO i = 1, kdlon
463	PTL(i,k) = (t(iof+i,k)+t(iof+i,k-1))*0.5
464	ENDDO
465	ENDDO
466	DO k = 1, kflev
467	DO i = 1, kdlon
468	PDP(i,k) = paprs(iof+i,k)-paprs(iof+i,k+1)
469	PTAVE(i,k) = t(iof+i,k)
470	PWV(i,k) = MAX (q(iof+i,k), 1.0e-12)
471	PQS(i,k) = PWV(i,k)
472	! Confert from column density of ozone in a cell, in kDU, to a mass fraction
473	POZON(i,k, :) = wo(iof+i, k, :) * RG * dobson_u * 1e3 &
474	/ (paprs(iof+i, k) - paprs(iof+i, k+1))
475	! A activer pour CCMVAL on prend l'ozone impose (MPL 07042010)
476	! POZON(i,k,:) = wo(i,k,:)
477	! print *,'RADLWSW: POZON',k, POZON(i,k,1)
478	PCLDLD(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)
479	PCLDLU(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)
480	PCLDSW(i,k) = cldfra(iof+i,k)
481	PTAU(i,1,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! 1e-12 serait instable
482	PTAU(i,2,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! pour 32-bit machines
483	POMEGA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAU(i,1,k))
484	POMEGA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAU(i,2,k))
485	PCG(i,1,k) = 0.865
486	PCG(i,2,k) = 0.910
487	!-
488	! Introduced for aerosol indirect forcings.
489	! The following values use the cloud optical thickness calculated from
490	! present-day aerosol concentrations whereas the quantities without the
491	! "A" at the end are for pre-industial (natural-only) aerosol concentrations
492	!
493	PTAUA(i,1,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! 1e-12 serait instable
494	PTAUA(i,2,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! pour 32-bit machines
495	POMEGAA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAUA(i,1,k))
496	POMEGAA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAUA(i,2,k))
497	ENDDO
498	ENDDO
499
500	IF (type_trac == 'repr') THEN
501	#ifdef REPROBUS
502	ndimozon = size(wo, 3)
503	CALL RAD_INTERACTIF(POZON,iof)
504	#endif
505	END IF
506
507	!
508	DO k = 1, kflev+1
509	DO i = 1, kdlon
510	PPMB(i,k) = paprs(iof+i,k)/100.0
511	ENDDO
512	ENDDO
513	!
514	!!!!! Modif MPL 6.01.09 avec RRTM, on passe de 5 a 6
515	DO kk = 1, 6
516	DO k = 1, kflev
517	DO i = 1, kdlon
518	PAER(i,k,kk) = 1.0E-15 !!!!! A REVOIR (MPL)
519	ENDDO
520	ENDDO
521	ENDDO
522	DO k = 1, kflev
523	DO i = 1, kdlon
524	tauaero(i,k,:,1)=tau_aero(iof+i,k,:,1)
525	pizaero(i,k,:,1)=piz_aero(iof+i,k,:,1)
526	cgaero(i,k,:,1) =cg_aero(iof+i,k,:,1)
527	tauaero(i,k,:,2)=tau_aero(iof+i,k,:,2)
528	pizaero(i,k,:,2)=piz_aero(iof+i,k,:,2)
529	cgaero(i,k,:,2) =cg_aero(iof+i,k,:,2)
530	ENDDO
531	ENDDO
532
533	!
534	!===== iflag_rrtm ================================================
535	!
536	IF (iflag_rrtm == 0) THEN !!!! remettre 0 juste pour tester l'ancien rayt via rrtm
537	!--- Mise a zero des tableaux output du rayonnement LW-AR4 ----------
538	DO k = 1, kflev+1
539	DO i = 1, kdlon
540	! print *,'RADLWSW: boucle mise a zero i k',i,k
541	ZFLUP(i,k)=0.
542	ZFLDN(i,k)=0.
543	ZFLUP0(i,k)=0.
544	ZFLDN0(i,k)=0.
545	ZLWFT0_i(i,k)=0.
546	ZFLUCUP_i(i,k)=0.
547	ZFLUCDWN_i(i,k)=0.
548	ENDDO
549	ENDDO
550	DO k = 1, kflev
551	DO i = 1, kdlon
552	zcool(i,k)=0.
553	zcool0(i,k)=0.
554	ENDDO
555	ENDDO
556	DO i = 1, kdlon
557	ztoplw(i)=0.
558	zsollw(i)=0.
559	ztoplw0(i)=0.
560	zsollw0(i)=0.
561	zsollwdown(i)=0.
562	ENDDO
563	! Old radiation scheme, used for AR4 runs
564	! average day-night ozone for longwave
565	CALL LW_LMDAR4(&
566	PPMB, PDP,&
567	PPSOL,PDT0,PEMIS,&
568	PTL, PTAVE, PWV, POZON(:, :, 1), PAER,&
569	PCLDLD,PCLDLU,&
570	PVIEW,&
571	zcool, zcool0,&
572	ztoplw,zsollw,ztoplw0,zsollw0,&
573	zsollwdown,&
574	ZFLUP, ZFLDN, ZFLUP0,ZFLDN0)
575	!----- Mise a zero des tableaux output du rayonnement SW-AR4
576	DO k = 1, kflev+1
577	DO i = 1, kdlon
578	ZFSUP(i,k)=0.
579	ZFSDN(i,k)=0.
580	ZFSUP0(i,k)=0.
581	ZFSDN0(i,k)=0.
582	ZSWFT0_i(i,k)=0.
583	ZFCUP_i(i,k)=0.
584	ZFCDWN_i(i,k)=0.
585	ENDDO
586	ENDDO
587	DO k = 1, kflev
588	DO i = 1, kdlon
589	zheat(i,k)=0.
590	zheat0(i,k)=0.
591	ENDDO
592	ENDDO
593	DO i = 1, kdlon
594	zalbpla(i)=0.
595	ztopsw(i)=0.
596	zsolsw(i)=0.
597	ztopsw0(i)=0.
598	zsolsw0(i)=0.
599	ztopswadaero(i)=0.
600	zsolswadaero(i)=0.
601	ztopswaiaero(i)=0.
602	zsolswaiaero(i)=0.
603	ENDDO
604	! print *,'Avant SW_LMDAR4: PSCT zrmu0 zfract',PSCT, zrmu0, zfract
605	! daylight ozone, if we have it, for short wave
606	IF (.NOT. new_aod) THEN
607	! use old version
608	CALL SW_LMDAR4(PSCT, zrmu0, zfract,&
609	PPMB, PDP, &
610	PPSOL, PALBD, PALBP,&
611	PTAVE, PWV, PQS, POZON(:, :, size(wo, 3)), PAER,&
612	PCLDSW, PTAU, POMEGA, PCG,&
613	zheat, zheat0,&
614	zalbpla,ztopsw,zsolsw,ztopsw0,zsolsw0,&
615	ZFSUP,ZFSDN,ZFSUP0,ZFSDN0,&
616	tauaero(:,:,5,:), pizaero(:,:,5,:), cgaero(:,:,5,:),&
617	PTAUA, POMEGAA,&
618	ztopswadaero,zsolswadaero,&
619	ztopswaiaero,zsolswaiaero,&
620	ok_ade, ok_aie)
621
622	ELSE ! new_aod=T
623	CALL SW_AEROAR4(PSCT, zrmu0, zfract,&
624	PPMB, PDP,&
625	PPSOL, PALBD, PALBP,&
626	PTAVE, PWV, PQS, POZON(:, :, size(wo, 3)), PAER,&
627	PCLDSW, PTAU, POMEGA, PCG,&
628	zheat, zheat0,&
629	zalbpla,ztopsw,zsolsw,ztopsw0,zsolsw0,&
630	ZFSUP,ZFSDN,ZFSUP0,ZFSDN0,&
631	tauaero, pizaero, cgaero, &
632	PTAUA, POMEGAA,&
633	ztopswadaero,zsolswadaero,&
634	ztopswad0aero,zsolswad0aero,&
635	ztopswaiaero,zsolswaiaero, &
636	ztopsw_aero,ztopsw0_aero,&
637	zsolsw_aero,zsolsw0_aero,&
638	ztopswcf_aero,zsolswcf_aero, &
639	ok_ade, ok_aie, flag_aerosol,flag_aerosol_strat)
640	ENDIF
641
642	ZSWFT0_i(:,:) = ZFSDN0(:,:)-ZFSUP0(:,:)
643	ZLWFT0_i(:,:) =-ZFLDN0(:,:)-ZFLUP0(:,:)
644
645	DO i=1,kdlon
646	DO k=1,kflev+1
647	! print *,'iof i k klon klev=',iof,i,k,klon,klev
648	lwdn0 ( iof+i,k) = ZFLDN0 ( i,k)
649	lwdn ( iof+i,k) = ZFLDN ( i,k)
650	lwup0 ( iof+i,k) = ZFLUP0 ( i,k)
651	lwup ( iof+i,k) = ZFLUP ( i,k)
652	swdn0 ( iof+i,k) = ZFSDN0 ( i,k)
653	swdn ( iof+i,k) = ZFSDN ( i,k)
654	swup0 ( iof+i,k) = ZFSUP0 ( i,k)
655	swup ( iof+i,k) = ZFSUP ( i,k)
656	ENDDO
657	ENDDO
658	! print*,'SW_AR4 ZFSDN0 1 , klev:',ZFSDN0(1:klon,1),ZFSDN0(1:klon,klev)
659	! print*,'SW_AR4 swdn0 1 , klev:',swdn0(1:klon,1),swdn0(1:klon,klev)
660	! print*,'SW_AR4 ZFSUP0 1 , klev:',ZFSUP0(1:klon,1),ZFSUP0(1:klon,klev)
661	! print*,'SW_AR4 swup0 1 , klev:',swup0(1:klon,1),swup0(1:klon,klev)
662	! print*,'SW_AR4 ZFSDN 1 , klev:',ZFSDN(1:klon,1) ,ZFSDN(1:klon,klev)
663	! print*,'SW_AR4 ZFSUP 1 , klev:',ZFSUP(1:klon,1) ,ZFSUP(1:klon,klev)
664	ELSE
665	#ifdef CPP_RRTM
666	! if (prt_level.gt.10)write(lunout,*)'CPP_RRTM=.T.'
667	!===== iflag_rrtm=1, on passe dans SW via RECMWFL ===============
668
669	DO k = 1, kflev+1
670	DO i = 1, kdlon
671	ZEMTD_i(i,k)=0.
672	ZEMTU_i(i,k)=0.
673	ZTRSO_i(i,k)=0.
674	ZTH_i(i,k)=0.
675	ZLWFT_i(i,k)=0.
676	ZSWFT_i(i,k)=0.
677	ZFLUX_i(i,1,k)=0.
678	ZFLUX_i(i,2,k)=0.
679	ZFLUC_i(i,1,k)=0.
680	ZFLUC_i(i,2,k)=0.
681	ZFSDWN_i(i,k)=0.
682	ZFCDWN_i(i,k)=0.
683	ZFSUP_i(i,k)=0.
684	ZFCUP_i(i,k)=0.
685	ENDDO
686	ENDDO
687	!
688	!--OB
689	!--aerosol TOT - anthropogenic+natural
690	!--aerosol NAT - natural only
691	!
692	DO i = 1, kdlon
693	DO k = 1, kflev
694	DO kk=1, NSW
695	!
696	PTAU_TOT(i,kflev+1-k,kk)=tau_aero_sw_rrtm(i,k,2,kk)
697	PPIZA_TOT(i,kflev+1-k,kk)=piz_aero_sw_rrtm(i,k,2,kk)
698	PCGA_TOT(i,kflev+1-k,kk)=cg_aero_sw_rrtm(i,k,2,kk)
699	!
700	PTAU_NAT(i,kflev+1-k,kk)=tau_aero_sw_rrtm(i,k,1,kk)
701	PPIZA_NAT(i,kflev+1-k,kk)=piz_aero_sw_rrtm(i,k,1,kk)
702	PCGA_NAT(i,kflev+1-k,kk)=cg_aero_sw_rrtm(i,k,1,kk)
703	!
704	ENDDO
705	ENDDO
706	ENDDO
707	!-end OB
708	!
709	!--C. Kleinschmitt
710	!--aerosol TOT - anthropogenic+natural
711	!--aerosol NAT - natural only
712	!
713	DO i = 1, kdlon
714	DO k = 1, kflev
715	DO kk=1, NLW
716	!
717	PTAU_LW_TOT(i,kflev+1-k,kk)=tau_aero_lw_rrtm(i,k,2,kk)
718	PTAU_LW_NAT(i,kflev+1-k,kk)=tau_aero_lw_rrtm(i,k,1,kk)
719	!
720	ENDDO
721	ENDDO
722	ENDDO
723	!-end C. Kleinschmitt
724	!
725	DO i = 1, kdlon
726	ZCTRSO(i,1)=0.
727	ZCTRSO(i,2)=0.
728	ZCEMTR(i,1)=0.
729	ZCEMTR(i,2)=0.
730	ZTRSOD(i)=0.
731	ZLWFC(i,1)=0.
732	ZLWFC(i,2)=0.
733	ZSWFC(i,1)=0.
734	ZSWFC(i,2)=0.
735	PFSDNN(i)=0.
736	PFSDNV(i)=0.
737	DO kk = 1, NSW
738	PSFSWDIR(i,kk)=0.
739	PSFSWDIF(i,kk)=0.
740	ENDDO
741	ENDDO
742	!----- Fin des mises a zero des tableaux output de RECMWF -------------------
743	! GEMU(1:klon)=sin(rlatd(1:klon))
744	! On met les donnees dans l'ordre des niveaux arpege
745	paprs_i(:,1)=paprs(:,klev+1)
746	do k=1,klev
747	paprs_i(1:klon,k+1) =paprs(1:klon,klev+1-k)
748	pplay_i(1:klon,k) =pplay(1:klon,klev+1-k)
749	cldfra_i(1:klon,k) =cldfra(1:klon,klev+1-k)
750	PDP_i(1:klon,k) =PDP(1:klon,klev+1-k)
751	t_i(1:klon,k) =t(1:klon,klev+1-k)
752	q_i(1:klon,k) =q(1:klon,klev+1-k)
753	qsat_i(1:klon,k) =qsat(1:klon,klev+1-k)
754	flwc_i(1:klon,k) =flwc(1:klon,klev+1-k)
755	fiwc_i(1:klon,k) =fiwc(1:klon,klev+1-k)
756	ref_liq_i(1:klon,k) =ref_liq(1:klon,klev+1-k)
757	ref_ice_i(1:klon,k) =ref_ice(1:klon,klev+1-k)
758	!-OB
759	ref_liq_pi_i(1:klon,k) =ref_liq_pi(1:klon,klev+1-k)
760	ref_ice_pi_i(1:klon,k) =ref_ice_pi(1:klon,klev+1-k)
761	enddo
762	do k=1,kflev
763	POZON_i(1:klon,k,:)=POZON(1:klon,kflev+1-k,:)
764	!!! POZON_i(1:klon,k)=POZON(1:klon,k) !!! on laisse 1=sol et klev=top
765	! print *,'Juste avant RECMWFL: k tsol temp',k,tsol,t(1,k)
766	!!!!!!! Modif MPL 6.01.09 avec RRTM, on passe de 5 a 6
767	do i=1,6
768	PAER_i(1:klon,k,i)=PAER(1:klon,kflev+1-k,i)
769	enddo
770	enddo
771	! print *,'RADLWSW: avant RECMWFL, RI0,rmu0=',solaire,rmu0
772
773	! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
774	! La version ARPEGE1D utilise differentes valeurs de la constante
775	! solaire suivant le rayonnement utilise.
776	! A controler ...
777	! SOLAR FLUX AT THE TOP (/YOMPHY3/)
778	! introduce season correction
779	!--------------------------------------
780	! RII0 = RIP0
781	! IF(LRAYFM)
782	! RII0 = RIP0M ! =rip0m if Morcrette non-each time step call.
783	! IF(LRAYFM15)
784	! RII0 = RIP0M15 ! =rip0m if Morcrette non-each time step call.
785	RII0=solaire/zdist/zdist
786	!print*,'+++ radlwsw: solaire ,RII0',solaire,RII0
787	! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
788	! Ancien appel a RECMWF (celui du cy25)
789	! CALL RECMWF (ist , iend, klon , ktdia , klev , kmode ,
790	! s PALBD , PALBP , paprs_i , pplay_i , RCO2 , cldfra_i,
791	! s POZON_i , PAER_i , PDP_i , PEMIS , GEMU , rmu0,
792	! s q_i , qsat_i , fiwc_i , flwc_i , zmasq , t_i ,tsol,
793	! s ZEMTD_i , ZEMTU_i , ZTRSO_i ,
794	! s ZTH_i , ZCTRSO , ZCEMTR , ZTRSOD ,
795	! s ZLWFC , ZLWFT_i , ZSWFC , ZSWFT_i ,
796	! s ZFLUX_i , ZFLUC_i , ZFSDWN_i, ZFSUP_i , ZFCDWN_i,ZFCUP_i)
797	! s 'RECMWF ')
798	!
799	if(lldebug) then
800	CALL writefield_phy('paprs_i',paprs_i,klev+1)
801	CALL writefield_phy('pplay_i',pplay_i,klev)
802	CALL writefield_phy('cldfra_i',cldfra_i,klev)
803	CALL writefield_phy('pozon_i',POZON_i,klev)
804	CALL writefield_phy('paer_i',PAER_i,klev)
805	CALL writefield_phy('pdp_i',PDP_i,klev)
806	CALL writefield_phy('q_i',q_i,klev)
807	CALL writefield_phy('qsat_i',qsat_i,klev)
808	CALL writefield_phy('fiwc_i',fiwc_i,klev)
809	CALL writefield_phy('flwc_i',flwc_i,klev)
810	CALL writefield_phy('t_i',t_i,klev)
811	CALL writefield_phy('palbd_new',PALBD_NEW,NSW)
812	CALL writefield_phy('palbp_new',PALBP_NEW,NSW)
813	endif
814
815	! Nouvel appel a RECMWF (celui du cy32t0)
816	CALL RECMWF_AERO (ist , iend, klon , ktdia , klev , kmode ,&
817	PALBD_NEW,PALBP_NEW, paprs_i , pplay_i , RCO2 , cldfra_i,&
818	POZON_i , PAER_i , PDP_i , PEMIS , rmu0 ,&
819	q_i , qsat_i , fiwc_i , flwc_i , zmasq , t_i ,tsol,&
820	ref_liq_i, ref_ice_i, &
821	ref_liq_pi_i, ref_ice_pi_i, & ! rajoute par OB pour diagnostiquer effet indirect
822	ZEMTD_i , ZEMTU_i , ZTRSO_i ,&
823	ZTH_i , ZCTRSO , ZCEMTR , ZTRSOD ,&
824	ZLWFC , ZLWFT_i , ZSWFC , ZSWFT_i ,&
825	PSFSWDIR , PSFSWDIF, PFSDNN , PFSDNV ,&
826	PPIZA_TOT, PCGA_TOT,PTAU_TOT,&
827	PPIZA_NAT, PCGA_NAT,PTAU_NAT, & ! rajoute par OB pour diagnostiquer effet direct
828	PTAU_LW_TOT, PTAU_LW_NAT, & ! rajoute par C. Kleinschmitt
829	ZFLUX_i , ZFLUC_i ,&
830	ZFSDWN_i , ZFSUP_i , ZFCDWN_i, ZFCUP_i,&
831	ZTOPSWADAERO,ZSOLSWADAERO,& ! rajoute par OB pour diagnostics
832	ZTOPSWAD0AERO,ZSOLSWAD0AERO,&
833	ZTOPSWAIAERO,ZSOLSWAIAERO, &
834	ZTOPSWCF_AERO,ZSOLSWCF_AERO, &
835	ZTOPLWADAERO,ZSOLLWADAERO,& ! rajoute par C. Kleinscmitt pour LW diagnostics
836	ZTOPLWAD0AERO,ZSOLLWAD0AERO,&
837	ZTOPLWAIAERO,ZSOLLWAIAERO, &
838	ok_ade, ok_aie, flag_aerosol,flag_aerosol_strat) ! flags aerosols
839
840	! print *,'RADLWSW: apres RECMWF'
841	if(lldebug) then
842	CALL writefield_phy('zemtd_i',ZEMTD_i,klev+1)
843	CALL writefield_phy('zemtu_i',ZEMTU_i,klev+1)
844	CALL writefield_phy('ztrso_i',ZTRSO_i,klev+1)
845	CALL writefield_phy('zth_i',ZTH_i,klev+1)
846	CALL writefield_phy('zctrso',ZCTRSO,2)
847	CALL writefield_phy('zcemtr',ZCEMTR,2)
848	CALL writefield_phy('ztrsod',ZTRSOD,1)
849	CALL writefield_phy('zlwfc',ZLWFC,2)
850	CALL writefield_phy('zlwft_i',ZLWFT_i,klev+1)
851	CALL writefield_phy('zswfc',ZSWFC,2)
852	CALL writefield_phy('zswft_i',ZSWFT_i,klev+1)
853	CALL writefield_phy('psfswdir',PSFSWDIR,6)
854	CALL writefield_phy('psfswdif',PSFSWDIF,6)
855	CALL writefield_phy('pfsdnn',PFSDNN,1)
856	CALL writefield_phy('pfsdnv',PFSDNV,1)
857	CALL writefield_phy('ppiza_dst',PPIZA_TOT,klev)
858	CALL writefield_phy('pcga_dst',PCGA_TOT,klev)
859	CALL writefield_phy('ptaurel_dst',PTAU_TOT,klev)
860	CALL writefield_phy('zflux_i',ZFLUX_i,klev+1)
861	CALL writefield_phy('zfluc_i',ZFLUC_i,klev+1)
862	CALL writefield_phy('zfsdwn_i',ZFSDWN_i,klev+1)
863	CALL writefield_phy('zfsup_i',ZFSUP_i,klev+1)
864	CALL writefield_phy('zfcdwn_i',ZFCDWN_i,klev+1)
865	CALL writefield_phy('zfcup_i',ZFCUP_i,klev+1)
866	endif
867	! --------- output RECMWFL
868	! ZEMTD (KPROMA,KLEV+1) ; TOTAL DOWNWARD LONGWAVE EMISSIVITY
869	! ZEMTU (KPROMA,KLEV+1) ; TOTAL UPWARD LONGWAVE EMISSIVITY
870	! ZTRSO (KPROMA,KLEV+1) ; TOTAL SHORTWAVE TRANSMISSIVITY
871	! ZTH (KPROMA,KLEV+1) ; HALF LEVEL TEMPERATURE
872	! ZCTRSO (KPROMA,2) ; CLEAR-SKY SHORTWAVE TRANSMISSIVITY
873	! ZCEMTR (KPROMA,2) ; CLEAR-SKY NET LONGWAVE EMISSIVITY
874	! ZTRSOD (KPROMA) ; TOTAL-SKY SURFACE SW TRANSMISSITY
875	! ZLWFC (KPROMA,2) ; CLEAR-SKY LONGWAVE FLUXES
876	! ZLWFT (KPROMA,KLEV+1) ; TOTAL-SKY LONGWAVE FLUXES
877	! ZSWFC (KPROMA,2) ; CLEAR-SKY SHORTWAVE FLUXES
878	! ZSWFT (KPROMA,KLEV+1) ; TOTAL-SKY SHORTWAVE FLUXES
879	! PPIZA_TOT (KPROMA,KLEV,NSW); Single scattering albedo of total aerosols
880	! PCGA_TOT (KPROMA,KLEV,NSW); Assymetry factor for total aerosols
881	! PTAU_TOT (KPROMA,KLEV,NSW); Optical depth of total aerosols
882	! PPIZA_NAT (KPROMA,KLEV,NSW); Single scattering albedo of natural aerosols
883	! PCGA_NAT (KPROMA,KLEV,NSW); Assymetry factor for natural aerosols
884	! PTAU_NAT (KPROMA,KLEV,NSW); Optical depth of natiral aerosols
885	! PTAU_LW_TOT (KPROMA,KLEV,NLW); LW Optical depth of total aerosols
886	! PTAU_LW_NAT (KPROMA,KLEV,NLW); LW Optical depth of natural aerosols
887	! PSFSWDIR (KPROMA,NSW) ;
888	! PSFSWDIF (KPROMA,NSW) ;
889	! PFSDNN (KPROMA) ;
890	! PFSDNV (KPROMA) ;
891	! ---------
892	! ---------
893	! On retablit l'ordre des niveaux lmd pour les tableaux de sortie
894	! D autre part, on multiplie les resultats SW par fract pour etre coherent
895	! avec l ancien rayonnement AR4. Si nuit, fract=0 donc pas de
896	! rayonnement SW. (MPL 260609)
897	DO k=0,klev
898	DO i=1,klon
899	ZEMTD(i,k+1) = ZEMTD_i(i,k+1)
900	ZEMTU(i,k+1) = ZEMTU_i(i,k+1)
901	ZTRSO(i,k+1) = ZTRSO_i(i,k+1)
902	ZTH(i,k+1) = ZTH_i(i,k+1)
903	! ZLWFT(i,k+1) = ZLWFT_i(i,klev+1-k)
904	! ZSWFT(i,k+1) = ZSWFT_i(i,klev+1-k)
905	ZFLUP(i,k+1) = ZFLUX_i(i,1,k+1)
906	ZFLDN(i,k+1) = ZFLUX_i(i,2,k+1)
907	ZFLUP0(i,k+1) = ZFLUC_i(i,1,k+1)
908	ZFLDN0(i,k+1) = ZFLUC_i(i,2,k+1)
909	ZFSDN(i,k+1) = ZFSDWN_i(i,k+1)*fract(i)
910	ZFSDN0(i,k+1) = ZFCDWN_i(i,k+1)*fract(i)
911	ZFSUP (i,k+1) = ZFSUP_i(i,k+1)*fract(i)
912	ZFSUP0(i,k+1) = ZFCUP_i(i,k+1)*fract(i)
913	! Nouveau calcul car visiblement ZSWFT et ZSWFC sont nuls dans RRTM cy32
914	! en sortie de radlsw.F90 - MPL 7.01.09
915	ZSWFT(i,k+1) = (ZFSDWN_i(i,k+1)-ZFSUP_i(i,k+1))*fract(i)
916	ZSWFT0_i(i,k+1) = (ZFCDWN_i(i,k+1)-ZFCUP_i(i,k+1))*fract(i)
917	! WRITE(*,'("FSDN FSUP FCDN FCUP: ",4E12.5)') ZFSDWN_i(i,k+1),&
918	! ZFSUP_i(i,k+1),ZFCDWN_i(i,k+1),ZFCUP_i(i,k+1)
919	ZLWFT(i,k+1) =-ZFLUX_i(i,2,k+1)-ZFLUX_i(i,1,k+1)
920	ZLWFT0_i(i,k+1)=-ZFLUC_i(i,2,k+1)-ZFLUC_i(i,1,k+1)
921	! print *,'FLUX2 FLUX1 FLUC2 FLUC1',ZFLUX_i(i,2,k+1),&
922	! & ZFLUX_i(i,1,k+1),ZFLUC_i(i,2,k+1),ZFLUC_i(i,1,k+1)
923	ENDDO
924	ENDDO
925
926	!--ajout OB
927	ZTOPSWADAERO(:) =ZTOPSWADAERO(:) *fract(:)
928	ZSOLSWADAERO(:) =ZSOLSWADAERO(:) *fract(:)
929	ZTOPSWAD0AERO(:)=ZTOPSWAD0AERO(:)*fract(:)
930	ZSOLSWAD0AERO(:)=ZSOLSWAD0AERO(:)*fract(:)
931	ZTOPSWAIAERO(:) =ZTOPSWAIAERO(:) *fract(:)
932	ZSOLSWAIAERO(:) =ZSOLSWAIAERO(:) *fract(:)
933	ZTOPSWCF_AERO(:,1)=ZTOPSWCF_AERO(:,1)*fract(:)
934	ZTOPSWCF_AERO(:,2)=ZTOPSWCF_AERO(:,2)*fract(:)
935	ZTOPSWCF_AERO(:,3)=ZTOPSWCF_AERO(:,3)*fract(:)
936	ZSOLSWCF_AERO(:,1)=ZSOLSWCF_AERO(:,1)*fract(:)
937	ZSOLSWCF_AERO(:,2)=ZSOLSWCF_AERO(:,2)*fract(:)
938	ZSOLSWCF_AERO(:,3)=ZSOLSWCF_AERO(:,3)*fract(:)
939
940	! print*,'SW_RRTM ZFSDN0 1 , klev:',ZFSDN0(1:klon,1),ZFSDN0(1:klon,klev)
941	! print*,'SW_RRTM ZFSUP0 1 , klev:',ZFSUP0(1:klon,1),ZFSUP0(1:klon,klev)
942	! print*,'SW_RRTM ZFSDN 1 , klev:',ZFSDN(1:klon,1),ZFSDN(1:klon,klev)
943	! print*,'SW_RRTM ZFSUP 1 , klev:',ZFSUP(1:klon,1),ZFSUP(1:klon,klev)
944	! print*,'OK1'
945	! ---------
946	! ---------
947	! On renseigne les champs LMDz, pour avoir la meme chose qu'en sortie de
948	! LW_LMDAR4 et SW_LMDAR4
949	DO i = 1, kdlon
950	zsolsw(i) = ZSWFT(i,1)
951	zsolsw0(i) = ZSWFT0_i(i,1)
952	! zsolsw0(i) = ZFSDN0(i,1) -ZFSUP0(i,1)
953	ztopsw(i) = ZSWFT(i,klev+1)
954	ztopsw0(i) = ZSWFT0_i(i,klev+1)
955	! ztopsw0(i) = ZFSDN0(i,klev+1)-ZFSUP0(i,klev+1)
956	!
957	! zsollw(i) = ZFLDN(i,1) -ZFLUP(i,1)
958	! zsollw0(i) = ZFLDN0(i,1) -ZFLUP0(i,1)
959	! ztoplw(i) = ZFLDN(i,klev+1) -ZFLUP(i,klev+1)
960	! ztoplw0(i) = ZFLDN0(i,klev+1)-ZFLUP0(i,klev+1)
961	zsollw(i) = ZLWFT(i,1)
962	zsollw0(i) = ZLWFT0_i(i,1)
963	ztoplw(i) = ZLWFT(i,klev+1)*(-1)
964	ztoplw0(i) = ZLWFT0_i(i,klev+1)*(-1)
965	!
966	IF (fract(i) == 0.) THEN
967	!!!!! A REVOIR MPL (20090630) ca n a pas de sens quand fract=0
968	! pas plus que dans le sw_AR4
969	zalbpla(i) = 1.0e+39
970	ELSE
971	zalbpla(i) = ZFSUP(i,klev+1)/ZFSDN(i,klev+1)
972	ENDIF
973	!!! 5 juin 2015
974	!!! Correction MP bug RRTM
975	zsollwdown(i)= -1.*ZFLDN(i,1)
976	ENDDO
977	! print*,'OK2'
978
979	! extrait de SW_AR4
980	! DO k = 1, KFLEV
981	! kpl1 = k+1
982	! DO i = 1, KDLON
983	! PHEAT(i,k) = -(ZFSUP(i,kpl1)-ZFSUP(i,k)) -(ZFSDN(i,k)-ZFSDN(i,kpl1))
984	! PHEAT(i,k) = PHEAT(i,k) * RDAY*RG/RCPD / PDP(i,k)
985	! ZLWFT(klon,k),ZSWFT
986
987	do k=1,kflev
988	do i=1,kdlon
989	zheat(i,k)=(ZSWFT(i,k+1)-ZSWFT(i,k))RDAYRG/RCPD/PDP(i,k)
990	zheat0(i,k)=(ZSWFT0_i(i,k+1)-ZSWFT0_i(i,k))RDAYRG/RCPD/PDP(i,k)
991	zcool(i,k)=(ZLWFT(i,k)-ZLWFT(i,k+1))RDAYRG/RCPD/PDP(i,k)
992	zcool0(i,k)=(ZLWFT0_i(i,k)-ZLWFT0_i(i,k+1))RDAYRG/RCPD/PDP(i,k)
993	! print *,'heat cool heat0 cool0 ',zheat(i,k),zcool(i,k),zheat0(i,k),zcool0(i,k)
994	! ZFLUCUP_i(i,k)=ZFLUC_i(i,1,k)
995	! ZFLUCDWN_i(i,k)=ZFLUC_i(i,2,k)
996	enddo
997	enddo
998	#else
999	abort_message="You should compile with -rrtm if running with iflag_rrtm=1"
1000	call abort_physic(modname, abort_message, 1)
1001	#endif
1002	ENDIF ! iflag_rrtm
1003	!======================================================================
1004
1005	DO i = 1, kdlon
1006	topsw(iof+i) = ztopsw(i)
1007	toplw(iof+i) = ztoplw(i)
1008	solsw(iof+i) = zsolsw(i)
1009	sollw(iof+i) = zsollw(i)
1010	sollwdown(iof+i) = zsollwdown(i)
1011	DO k = 1, kflev+1
1012	lwdn0 ( iof+i,k) = ZFLDN0 ( i,k)
1013	lwdn ( iof+i,k) = ZFLDN ( i,k)
1014	lwup0 ( iof+i,k) = ZFLUP0 ( i,k)
1015	lwup ( iof+i,k) = ZFLUP ( i,k)
1016	ENDDO
1017	topsw0(iof+i) = ztopsw0(i)
1018	toplw0(iof+i) = ztoplw0(i)
1019	solsw0(iof+i) = zsolsw0(i)
1020	sollw0(iof+i) = zsollw0(i)
1021	albpla(iof+i) = zalbpla(i)
1022
1023	DO k = 1, kflev+1
1024	swdn0 ( iof+i,k) = ZFSDN0 ( i,k)
1025	swdn ( iof+i,k) = ZFSDN ( i,k)
1026	swup0 ( iof+i,k) = ZFSUP0 ( i,k)
1027	swup ( iof+i,k) = ZFSUP ( i,k)
1028	ENDDO
1029	ENDDO
1030	!-transform the aerosol forcings, if they have
1031	! to be calculated
1032	IF (ok_ade) THEN
1033	DO i = 1, kdlon
1034	topswad_aero(iof+i) = ztopswadaero(i)
1035	topswad0_aero(iof+i) = ztopswad0aero(i)
1036	solswad_aero(iof+i) = zsolswadaero(i)
1037	solswad0_aero(iof+i) = zsolswad0aero(i)
1038	! MS the following lines seem to be wrong, why is iof on right hand side???
1039	! topsw_aero(iof+i,:) = ztopsw_aero(iof+i,:)
1040	! topsw0_aero(iof+i,:) = ztopsw0_aero(iof+i,:)
1041	! solsw_aero(iof+i,:) = zsolsw_aero(iof+i,:)
1042	! solsw0_aero(iof+i,:) = zsolsw0_aero(iof+i,:)
1043	topsw_aero(iof+i,:) = ztopsw_aero(i,:)
1044	topsw0_aero(iof+i,:) = ztopsw0_aero(i,:)
1045	solsw_aero(iof+i,:) = zsolsw_aero(i,:)
1046	solsw0_aero(iof+i,:) = zsolsw0_aero(i,:)
1047	topswcf_aero(iof+i,:) = ztopswcf_aero(i,:)
1048	solswcf_aero(iof+i,:) = zsolswcf_aero(i,:)
1049	!-LW
1050	toplwad_aero(iof+i) = ztoplwadaero(i)
1051	toplwad0_aero(iof+i) = ztoplwad0aero(i)
1052	sollwad_aero(iof+i) = zsollwadaero(i)
1053	sollwad0_aero(iof+i) = zsollwad0aero(i)
1054	ENDDO
1055	ELSE
1056	DO i = 1, kdlon
1057	topswad_aero(iof+i) = 0.0
1058	solswad_aero(iof+i) = 0.0
1059	topswad0_aero(iof+i) = 0.0
1060	solswad0_aero(iof+i) = 0.0
1061	topsw_aero(iof+i,:) = 0.
1062	topsw0_aero(iof+i,:) =0.
1063	solsw_aero(iof+i,:) = 0.
1064	solsw0_aero(iof+i,:) = 0.
1065	!-LW
1066	toplwad_aero(iof+i) = 0.0
1067	sollwad_aero(iof+i) = 0.0
1068	toplwad0_aero(iof+i) = 0.0
1069	sollwad0_aero(iof+i) = 0.0
1070	ENDDO
1071	ENDIF
1072	IF (ok_aie) THEN
1073	DO i = 1, kdlon
1074	topswai_aero(iof+i) = ztopswaiaero(i)
1075	solswai_aero(iof+i) = zsolswaiaero(i)
1076	!-LW
1077	toplwai_aero(iof+i) = ztoplwaiaero(i)
1078	sollwai_aero(iof+i) = zsollwaiaero(i)
1079	ENDDO
1080	ELSE
1081	DO i = 1, kdlon
1082	topswai_aero(iof+i) = 0.0
1083	solswai_aero(iof+i) = 0.0
1084	!-LW
1085	toplwai_aero(iof+i) = 0.0
1086	sollwai_aero(iof+i) = 0.0
1087	ENDDO
1088	ENDIF
1089	DO k = 1, kflev
1090	DO i = 1, kdlon
1091	! scale factor to take into account the difference between
1092	! dry air and watter vapour scpecifi! heat capacity
1093	zznormcp=1.0+RVTMP2*PWV(i,k)
1094	heat(iof+i,k) = zheat(i,k)/zznormcp
1095	cool(iof+i,k) = zcool(i,k)/zznormcp
1096	heat0(iof+i,k) = zheat0(i,k)/zznormcp
1097	cool0(iof+i,k) = zcool0(i,k)/zznormcp
1098	ENDDO
1099	ENDDO
1100
1101	ENDDO ! j = 1, nb_gr
1102
1103	END SUBROUTINE radlwsw
1104
1105	end module radlwsw_m

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

Download in other formats:

Original Format