Context Navigation

radlwsw_m.F90 @ 2390

Last change on this file since 2390 was 2366, checked in by oboucher, 10 years ago
Second batch of changes for diagnosing SW radiative diagnostics at the physics timestep resolution radsol is now computed in physiq.F90
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:keywords set to `Author Date Id Revi`
File size: 45.4 KB

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

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

Download in other formats:

Original Format