Context Navigation

radlwsw_m.F90 @ 2394

Last change on this file since 2394 was 2394, checked in by acozic, 9 years ago

to fit with inca we need to increase dimension of naero_grp in aero_mod.
When lmdz is running without inca dim 10 to 12 will be empty and there is no change in results

these modifications only impact the old physics coupled with inca

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 2394 2015-11-18 11:41:49Z acozic $
3	!
4	module radlwsw_m
5
6	IMPLICIT NONE
7
8	contains
9
10	SUBROUTINE radlwsw( &
11	dist, rmu0, fract, &
12	!albedo SB >>>
13	! paprs, pplay,tsol,alb1, alb2, &
14	paprs, pplay,tsol,SFRWL,alb_dir, alb_dif, &
15	!albedo SB <<<
16	t,q,wo,&
17	cldfra, cldemi, cldtaupd,&
18	ok_ade, ok_aie, 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 (T/F)
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	LOGICAL, 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	!
388	!-------------------------------------------
389	nb_gr = KLON / kdlon
390	IF (nb_gr*kdlon .NE. KLON) THEN
391	PRINT*, "kdlon mauvais:", KLON, kdlon, nb_gr
392	call abort_physic("radlwsw", "", 1)
393	ENDIF
394	IF (kflev .NE. KLEV) THEN
395	PRINT*, "kflev differe de KLEV, kflev, KLEV"
396	call abort_physic("radlwsw", "", 1)
397	ENDIF
398	!-------------------------------------------
399	DO k = 1, KLEV
400	DO i = 1, KLON
401	heat(i,k)=0.
402	cool(i,k)=0.
403	heat0(i,k)=0.
404	cool0(i,k)=0.
405	ENDDO
406	ENDDO
407	!
408	zdist = dist
409	!
410	PSCT = solaire/zdist/zdist
411
412	IF (type_trac == 'repr') THEN
413	#ifdef REPROBUS
414	if(ok_SUNTIME) PSCT = solaireTIME/zdist/zdist
415	print,'Constante solaire: ',PSCTzdist*zdist
416	#endif
417	END IF
418
419	DO j = 1, nb_gr
420	iof = kdlon*(j-1)
421	DO i = 1, kdlon
422	zfract(i) = fract(iof+i)
423	! zfract(i) = 1. !!!!!! essai MPL 19052010
424	zrmu0(i) = rmu0(iof+i)
425
426
427	!albedo SB >>>
428	! PALBD(i,1) = alb1(iof+i)
429	! PALBD(i,2) = alb2(iof+i)
430	! PALBD_NEW(i,1) = alb1(iof+i) !!!!! A REVOIR (MPL) PALBD_NEW en
431	! fonction bdes SW
432	! do kk=2,NSW
433	! PALBD_NEW(i,kk) = alb2(iof+i)
434	! enddo
435	! PALBP(i,1) = alb1(iof+i)
436	! PALBP(i,2) = alb2(iof+i)
437	!
438	! PALBP_NEW(i,1) = alb1(iof+i) !!!!! A REVOIR (MPL) PALBP_NEW en
439	! fonction bdes SW
440	! do kk=2,NSW
441	! PALBP_NEW(i,kk) = alb2(iof+i)
442	! enddo
443
444	if(iflag_rrtm==0)then
445	select case(nsw)
446	case(2)
447	PALBD(i,1)=alb_dif(iof+i,1)
448	PALBD(i,2)=alb_dif(iof+i,2)
449	PALBP(i,1)=alb_dir(iof+i,1)
450	PALBP(i,2)=alb_dir(iof+i,2)
451	case(4)
452	PALBD(i,1)=alb_dif(iof+i,1)
453	PALBD(i,2)=(alb_dif(iof+i,2)SFRWL(2)+alb_dif(iof+i,3)SFRWL(3) &
454	+alb_dif(iof+i,4)*SFRWL(4))/(SFRWL(2)+SFRWL(3)+SFRWL(4))
455	PALBP(i,1)=alb_dir(iof+i,1)
456	PALBP(i,2)=(alb_dir(iof+i,2)SFRWL(2)+alb_dir(iof+i,3)SFRWL(3) &
457	+alb_dir(iof+i,4)*SFRWL(4))/(SFRWL(2)+SFRWL(3)+SFRWL(4))
458	case(6)
459	PALBD(i,1)=(alb_dif(iof+i,1)SFRWL(1)+alb_dif(iof+i,2)SFRWL(2) &
460	+alb_dif(iof+i,3)*SFRWL(3))/(SFRWL(1)+SFRWL(2)+SFRWL(3))
461	PALBD(i,2)=(alb_dif(iof+i,4)SFRWL(4)+alb_dif(iof+i,5)SFRWL(5) &
462	+alb_dif(iof+i,6)*SFRWL(6))/(SFRWL(4)+SFRWL(5)+SFRWL(6))
463	PALBP(i,1)=(alb_dir(iof+i,1)SFRWL(1)+alb_dir(iof+i,2)SFRWL(2) &
464	+alb_dir(iof+i,3)*SFRWL(3))/(SFRWL(1)+SFRWL(2)+SFRWL(3))
465	PALBP(i,2)=(alb_dir(iof+i,4)SFRWL(4)+alb_dir(iof+i,5)SFRWL(5) &
466	+alb_dir(iof+i,6)*SFRWL(6))/(SFRWL(4)+SFRWL(5)+SFRWL(6))
467	end select
468	elseif(iflag_rrtm==1)then
469	DO kk=1,NSW
470	PALBD_NEW(i,kk)=alb_dif(iof+i,kk)
471	PALBP_NEW(i,kk)=alb_dir(iof+i,kk)
472	ENDDO
473	endif
474	!albedo SB <<<
475
476
477
478
479	PEMIS(i) = 1.0 !!!!! A REVOIR (MPL)
480	PVIEW(i) = 1.66
481	PPSOL(i) = paprs(iof+i,1)
482	zx_alpha1 = (paprs(iof+i,1)-pplay(iof+i,2))/(pplay(iof+i,1)-pplay(iof+i,2))
483	zx_alpha2 = 1.0 - zx_alpha1
484	PTL(i,1) = t(iof+i,1) * zx_alpha1 + t(iof+i,2) * zx_alpha2
485	PTL(i,KLEV+1) = t(iof+i,KLEV)
486	PDT0(i) = tsol(iof+i) - PTL(i,1)
487	ENDDO
488	DO k = 2, kflev
489	DO i = 1, kdlon
490	PTL(i,k) = (t(iof+i,k)+t(iof+i,k-1))*0.5
491	ENDDO
492	ENDDO
493	DO k = 1, kflev
494	DO i = 1, kdlon
495	PDP(i,k) = paprs(iof+i,k)-paprs(iof+i,k+1)
496	PTAVE(i,k) = t(iof+i,k)
497	PWV(i,k) = MAX (q(iof+i,k), 1.0e-12)
498	PQS(i,k) = PWV(i,k)
499	POZON(i,k, :) = wo(iof+i, k, :) * RG * dobson_u * 1e3 &
500	/ (paprs(iof+i, k) - paprs(iof+i, k+1))
501	! A activer pour CCMVAL on prend l'ozone impose (MPL 07042010)
502	! POZON(i,k,:) = wo(i,k,:)
503	! print *,'RADLWSW: POZON',k, POZON(i,k,1)
504	PCLDLD(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)
505	PCLDLU(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)
506	PCLDSW(i,k) = cldfra(iof+i,k)
507	PTAU(i,1,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! 1e-12 serait instable
508	PTAU(i,2,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! pour 32-bit machines
509	POMEGA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAU(i,1,k))
510	POMEGA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAU(i,2,k))
511	PCG(i,1,k) = 0.865
512	PCG(i,2,k) = 0.910
513	!-
514	! Introduced for aerosol indirect forcings.
515	! The following values use the cloud optical thickness calculated from
516	! present-day aerosol concentrations whereas the quantities without the
517	! "A" at the end are for pre-industial (natural-only) aerosol concentrations
518	!
519	PTAUA(i,1,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! 1e-12 serait instable
520	PTAUA(i,2,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! pour 32-bit machines
521	POMEGAA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAUA(i,1,k))
522	POMEGAA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAUA(i,2,k))
523	ENDDO
524	ENDDO
525
526	IF (type_trac == 'repr') THEN
527	#ifdef REPROBUS
528	ndimozon = size(wo, 3)
529	CALL RAD_INTERACTIF(POZON,iof)
530	#endif
531	END IF
532
533	!
534	DO k = 1, kflev+1
535	DO i = 1, kdlon
536	PPMB(i,k) = paprs(iof+i,k)/100.0
537	ENDDO
538	ENDDO
539	!
540	!!!!! Modif MPL 6.01.09 avec RRTM, on passe de 5 a 6
541	DO kk = 1, 6
542	DO k = 1, kflev
543	DO i = 1, kdlon
544	PAER(i,k,kk) = 1.0E-15 !!!!! A REVOIR (MPL)
545	ENDDO
546	ENDDO
547	ENDDO
548	DO k = 1, kflev
549	DO i = 1, kdlon
550	tauaero(i,k,:,1)=tau_aero(iof+i,k,:,1)
551	pizaero(i,k,:,1)=piz_aero(iof+i,k,:,1)
552	cgaero(i,k,:,1) =cg_aero(iof+i,k,:,1)
553	tauaero(i,k,:,2)=tau_aero(iof+i,k,:,2)
554	pizaero(i,k,:,2)=piz_aero(iof+i,k,:,2)
555	cgaero(i,k,:,2) =cg_aero(iof+i,k,:,2)
556	ENDDO
557	ENDDO
558
559	!
560	!===== iflag_rrtm ================================================
561	!
562	IF (iflag_rrtm == 0) THEN !!!! remettre 0 juste pour tester l'ancien rayt via rrtm
563	!--- Mise a zero des tableaux output du rayonnement LW-AR4 ----------
564	DO k = 1, kflev+1
565	DO i = 1, kdlon
566	! print *,'RADLWSW: boucle mise a zero i k',i,k
567	ZFLUP(i,k)=0.
568	ZFLDN(i,k)=0.
569	ZFLUP0(i,k)=0.
570	ZFLDN0(i,k)=0.
571	ZLWFT0_i(i,k)=0.
572	ZFLUCUP_i(i,k)=0.
573	ZFLUCDWN_i(i,k)=0.
574	ENDDO
575	ENDDO
576	DO k = 1, kflev
577	DO i = 1, kdlon
578	zcool(i,k)=0.
579	zcool0(i,k)=0.
580	ENDDO
581	ENDDO
582	DO i = 1, kdlon
583	ztoplw(i)=0.
584	zsollw(i)=0.
585	ztoplw0(i)=0.
586	zsollw0(i)=0.
587	zsollwdown(i)=0.
588	ENDDO
589	! Old radiation scheme, used for AR4 runs
590	! average day-night ozone for longwave
591	CALL LW_LMDAR4(&
592	PPMB, PDP,&
593	PPSOL,PDT0,PEMIS,&
594	PTL, PTAVE, PWV, POZON(:, :, 1), PAER,&
595	PCLDLD,PCLDLU,&
596	PVIEW,&
597	zcool, zcool0,&
598	ztoplw,zsollw,ztoplw0,zsollw0,&
599	zsollwdown,&
600	ZFLUP, ZFLDN, ZFLUP0,ZFLDN0)
601	!----- Mise a zero des tableaux output du rayonnement SW-AR4
602	DO k = 1, kflev+1
603	DO i = 1, kdlon
604	ZFSUP(i,k)=0.
605	ZFSDN(i,k)=0.
606	ZFSUP0(i,k)=0.
607	ZFSDN0(i,k)=0.
608	ZSWFT0_i(i,k)=0.
609	ZFCUP_i(i,k)=0.
610	ZFCDWN_i(i,k)=0.
611	ENDDO
612	ENDDO
613	DO k = 1, kflev
614	DO i = 1, kdlon
615	zheat(i,k)=0.
616	zheat0(i,k)=0.
617	ENDDO
618	ENDDO
619	DO i = 1, kdlon
620	zalbpla(i)=0.
621	ztopsw(i)=0.
622	zsolsw(i)=0.
623	ztopsw0(i)=0.
624	zsolsw0(i)=0.
625	ztopswadaero(i)=0.
626	zsolswadaero(i)=0.
627	ztopswaiaero(i)=0.
628	zsolswaiaero(i)=0.
629	ENDDO
630	! print *,'Avant SW_LMDAR4: PSCT zrmu0 zfract',PSCT, zrmu0, zfract
631	! daylight ozone, if we have it, for short wave
632	IF (.NOT. new_aod) THEN
633	! use old version
634	CALL SW_LMDAR4(PSCT, zrmu0, zfract,&
635	PPMB, PDP, &
636	PPSOL, PALBD, PALBP,&
637	PTAVE, PWV, PQS, POZON(:, :, size(wo, 3)), PAER,&
638	PCLDSW, PTAU, POMEGA, PCG,&
639	zheat, zheat0,&
640	zalbpla,ztopsw,zsolsw,ztopsw0,zsolsw0,&
641	ZFSUP,ZFSDN,ZFSUP0,ZFSDN0,&
642	tauaero(:,:,5,:), pizaero(:,:,5,:), cgaero(:,:,5,:),&
643	PTAUA, POMEGAA,&
644	ztopswadaero,zsolswadaero,&
645	ztopswaiaero,zsolswaiaero,&
646	ok_ade, ok_aie)
647
648	ELSE ! new_aod=T
649	CALL SW_AEROAR4(PSCT, zrmu0, zfract,&
650	PPMB, PDP,&
651	PPSOL, PALBD, PALBP,&
652	PTAVE, PWV, PQS, POZON(:, :, size(wo, 3)), PAER,&
653	PCLDSW, PTAU, POMEGA, PCG,&
654	zheat, zheat0,&
655	zalbpla,ztopsw,zsolsw,ztopsw0,zsolsw0,&
656	ZFSUP,ZFSDN,ZFSUP0,ZFSDN0,&
657	tauaero, pizaero, cgaero, &
658	PTAUA, POMEGAA,&
659	ztopswadaero,zsolswadaero,&
660	ztopswad0aero,zsolswad0aero,&
661	ztopswaiaero,zsolswaiaero, &
662	ztopsw_aero,ztopsw0_aero,&
663	zsolsw_aero,zsolsw0_aero,&
664	ztopswcf_aero,zsolswcf_aero, &
665	ok_ade, ok_aie, flag_aerosol,flag_aerosol_strat)
666	ENDIF
667
668
669	DO i=1,kdlon
670	DO k=1,kflev+1
671	ZSWFT0_i(1:klon,k) = ZFSDN0(1:klon,k)-ZFSUP0(1:klon,k)
672	ZLWFT0_i(1:klon,k)=-ZFLDN0(1:klon,k)-ZFLUP0(1:klon,k)
673	! print *,'iof i k klon klev=',iof,i,k,klon,klev
674	lwdn0 ( iof+i,k) = ZFLDN0 ( i,k)
675	lwdn ( iof+i,k) = ZFLDN ( i,k)
676	lwup0 ( iof+i,k) = ZFLUP0 ( i,k)
677	lwup ( iof+i,k) = ZFLUP ( i,k)
678	swdn0 ( iof+i,k) = ZFSDN0 ( i,k)
679	swdn ( iof+i,k) = ZFSDN ( i,k)
680	swup0 ( iof+i,k) = ZFSUP0 ( i,k)
681	swup ( iof+i,k) = ZFSUP ( i,k)
682	ENDDO
683	ENDDO
684	! print*,'SW_AR4 ZFSDN0 1 , klev:',ZFSDN0(1:klon,1),ZFSDN0(1:klon,klev)
685	! print*,'SW_AR4 swdn0 1 , klev:',swdn0(1:klon,1),swdn0(1:klon,klev)
686	! print*,'SW_AR4 ZFSUP0 1 , klev:',ZFSUP0(1:klon,1),ZFSUP0(1:klon,klev)
687	! print*,'SW_AR4 swup0 1 , klev:',swup0(1:klon,1),swup0(1:klon,klev)
688	! print*,'SW_AR4 ZFSDN 1 , klev:',ZFSDN(1:klon,1) ,ZFSDN(1:klon,klev)
689	! print*,'SW_AR4 ZFSUP 1 , klev:',ZFSUP(1:klon,1) ,ZFSUP(1:klon,klev)
690	ELSE
691	#ifdef CPP_RRTM
692	! if (prt_level.gt.10)write(lunout,*)'CPP_RRTM=.T.'
693	!===== iflag_rrtm=1, on passe dans SW via RECMWFL ===============
694
695	DO k = 1, kflev+1
696	DO i = 1, kdlon
697	ZEMTD_i(i,k)=0.
698	ZEMTU_i(i,k)=0.
699	ZTRSO_i(i,k)=0.
700	ZTH_i(i,k)=0.
701	ZLWFT_i(i,k)=0.
702	ZSWFT_i(i,k)=0.
703	ZFLUX_i(i,1,k)=0.
704	ZFLUX_i(i,2,k)=0.
705	ZFLUC_i(i,1,k)=0.
706	ZFLUC_i(i,2,k)=0.
707	ZFSDWN_i(i,k)=0.
708	ZFCDWN_i(i,k)=0.
709	ZFSUP_i(i,k)=0.
710	ZFCUP_i(i,k)=0.
711	ENDDO
712	ENDDO
713	!
714	!--OB
715	!--aerosol TOT - anthropogenic+natural
716	!--aerosol NAT - natural only
717	!
718	DO i = 1, kdlon
719	DO k = 1, kflev
720	DO kk=1, NSW
721	!
722	PTAU_TOT(i,kflev+1-k,kk)=tau_aero_sw_rrtm(i,k,2,kk)
723	PPIZA_TOT(i,kflev+1-k,kk)=piz_aero_sw_rrtm(i,k,2,kk)
724	PCGA_TOT(i,kflev+1-k,kk)=cg_aero_sw_rrtm(i,k,2,kk)
725	!
726	PTAU_NAT(i,kflev+1-k,kk)=tau_aero_sw_rrtm(i,k,1,kk)
727	PPIZA_NAT(i,kflev+1-k,kk)=piz_aero_sw_rrtm(i,k,1,kk)
728	PCGA_NAT(i,kflev+1-k,kk)=cg_aero_sw_rrtm(i,k,1,kk)
729	!
730	ENDDO
731	ENDDO
732	ENDDO
733	!-end OB
734	!
735	!--C. Kleinschmitt
736	!--aerosol TOT - anthropogenic+natural
737	!--aerosol NAT - natural only
738	!
739	DO i = 1, kdlon
740	DO k = 1, kflev
741	DO kk=1, NLW
742	!
743	PTAU_LW_TOT(i,kflev+1-k,kk)=tau_aero_lw_rrtm(i,k,2,kk)
744	PTAU_LW_NAT(i,kflev+1-k,kk)=tau_aero_lw_rrtm(i,k,1,kk)
745	!
746	ENDDO
747	ENDDO
748	ENDDO
749	!-end C. Kleinschmitt
750	!
751	DO i = 1, kdlon
752	ZCTRSO(i,1)=0.
753	ZCTRSO(i,2)=0.
754	ZCEMTR(i,1)=0.
755	ZCEMTR(i,2)=0.
756	ZTRSOD(i)=0.
757	ZLWFC(i,1)=0.
758	ZLWFC(i,2)=0.
759	ZSWFC(i,1)=0.
760	ZSWFC(i,2)=0.
761	PFSDNN(i)=0.
762	PFSDNV(i)=0.
763	DO kk = 1, NSW
764	PSFSWDIR(i,kk)=0.
765	PSFSWDIF(i,kk)=0.
766	ENDDO
767	ENDDO
768	!----- Fin des mises a zero des tableaux output de RECMWF -------------------
769	! GEMU(1:klon)=sin(rlatd(1:klon))
770	! On met les donnees dans l'ordre des niveaux arpege
771	paprs_i(:,1)=paprs(:,klev+1)
772	do k=1,klev
773	paprs_i(1:klon,k+1) =paprs(1:klon,klev+1-k)
774	pplay_i(1:klon,k) =pplay(1:klon,klev+1-k)
775	cldfra_i(1:klon,k) =cldfra(1:klon,klev+1-k)
776	PDP_i(1:klon,k) =PDP(1:klon,klev+1-k)
777	t_i(1:klon,k) =t(1:klon,klev+1-k)
778	q_i(1:klon,k) =q(1:klon,klev+1-k)
779	qsat_i(1:klon,k) =qsat(1:klon,klev+1-k)
780	flwc_i(1:klon,k) =flwc(1:klon,klev+1-k)
781	fiwc_i(1:klon,k) =fiwc(1:klon,klev+1-k)
782	ref_liq_i(1:klon,k) =ref_liq(1:klon,klev+1-k)
783	ref_ice_i(1:klon,k) =ref_ice(1:klon,klev+1-k)
784	!-OB
785	ref_liq_pi_i(1:klon,k) =ref_liq_pi(1:klon,klev+1-k)
786	ref_ice_pi_i(1:klon,k) =ref_ice_pi(1:klon,klev+1-k)
787	enddo
788	do k=1,kflev
789	POZON_i(1:klon,k,:)=POZON(1:klon,kflev+1-k,:)
790	!!! POZON_i(1:klon,k)=POZON(1:klon,k) !!! on laisse 1=sol et klev=top
791	! print *,'Juste avant RECMWFL: k tsol temp',k,tsol,t(1,k)
792	!!!!!!! Modif MPL 6.01.09 avec RRTM, on passe de 5 a 6
793	do i=1,6
794	PAER_i(1:klon,k,i)=PAER(1:klon,kflev+1-k,i)
795	enddo
796	enddo
797	! print *,'RADLWSW: avant RECMWFL, RI0,rmu0=',solaire,rmu0
798
799	! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
800	! La version ARPEGE1D utilise differentes valeurs de la constante
801	! solaire suivant le rayonnement utilise.
802	! A controler ...
803	! SOLAR FLUX AT THE TOP (/YOMPHY3/)
804	! introduce season correction
805	!--------------------------------------
806	! RII0 = RIP0
807	! IF(LRAYFM)
808	! RII0 = RIP0M ! =rip0m if Morcrette non-each time step call.
809	! IF(LRAYFM15)
810	! RII0 = RIP0M15 ! =rip0m if Morcrette non-each time step call.
811	RII0=solaire/zdist/zdist
812	!print*,'+++ radlwsw: solaire ,RII0',solaire,RII0
813	! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
814	! Ancien appel a RECMWF (celui du cy25)
815	! CALL RECMWF (ist , iend, klon , ktdia , klev , kmode ,
816	! s PALBD , PALBP , paprs_i , pplay_i , RCO2 , cldfra_i,
817	! s POZON_i , PAER_i , PDP_i , PEMIS , GEMU , rmu0,
818	! s q_i , qsat_i , fiwc_i , flwc_i , zmasq , t_i ,tsol,
819	! s ZEMTD_i , ZEMTU_i , ZTRSO_i ,
820	! s ZTH_i , ZCTRSO , ZCEMTR , ZTRSOD ,
821	! s ZLWFC , ZLWFT_i , ZSWFC , ZSWFT_i ,
822	! s ZFLUX_i , ZFLUC_i , ZFSDWN_i, ZFSUP_i , ZFCDWN_i,ZFCUP_i)
823	! s 'RECMWF ')
824	!
825	if(lldebug) then
826	CALL writefield_phy('paprs_i',paprs_i,klev+1)
827	CALL writefield_phy('pplay_i',pplay_i,klev)
828	CALL writefield_phy('cldfra_i',cldfra_i,klev)
829	CALL writefield_phy('pozon_i',POZON_i,klev)
830	CALL writefield_phy('paer_i',PAER_i,klev)
831	CALL writefield_phy('pdp_i',PDP_i,klev)
832	CALL writefield_phy('q_i',q_i,klev)
833	CALL writefield_phy('qsat_i',qsat_i,klev)
834	CALL writefield_phy('fiwc_i',fiwc_i,klev)
835	CALL writefield_phy('flwc_i',flwc_i,klev)
836	CALL writefield_phy('t_i',t_i,klev)
837	CALL writefield_phy('palbd_new',PALBD_NEW,NSW)
838	CALL writefield_phy('palbp_new',PALBP_NEW,NSW)
839	endif
840
841	! Nouvel appel a RECMWF (celui du cy32t0)
842	CALL RECMWF_AERO (ist , iend, klon , ktdia , klev , kmode ,&
843	PALBD_NEW,PALBP_NEW, paprs_i , pplay_i , RCO2 , cldfra_i,&
844	POZON_i , PAER_i , PDP_i , PEMIS , rmu0 ,&
845	q_i , qsat_i , fiwc_i , flwc_i , zmasq , t_i ,tsol,&
846	ref_liq_i, ref_ice_i, &
847	ref_liq_pi_i, ref_ice_pi_i, & ! rajoute par OB pour diagnostiquer effet indirect
848	ZEMTD_i , ZEMTU_i , ZTRSO_i ,&
849	ZTH_i , ZCTRSO , ZCEMTR , ZTRSOD ,&
850	ZLWFC , ZLWFT_i , ZSWFC , ZSWFT_i ,&
851	PSFSWDIR , PSFSWDIF, PFSDNN , PFSDNV ,&
852	PPIZA_TOT, PCGA_TOT,PTAU_TOT,&
853	PPIZA_NAT, PCGA_NAT,PTAU_NAT, & ! rajoute par OB pour diagnostiquer effet direct
854	PTAU_LW_TOT, PTAU_LW_NAT, & ! rajoute par C. Kleinschmitt
855	ZFLUX_i , ZFLUC_i ,&
856	ZFSDWN_i , ZFSUP_i , ZFCDWN_i, ZFCUP_i,&
857	ZTOPSWADAERO,ZSOLSWADAERO,& ! rajoute par OB pour diagnostics
858	ZTOPSWAD0AERO,ZSOLSWAD0AERO,&
859	ZTOPSWAIAERO,ZSOLSWAIAERO, &
860	ZTOPSWCF_AERO,ZSOLSWCF_AERO, &
861	ZTOPLWADAERO,ZSOLLWADAERO,& ! rajoute par C. Kleinscmitt pour LW diagnostics
862	ZTOPLWAD0AERO,ZSOLLWAD0AERO,&
863	ZTOPLWAIAERO,ZSOLLWAIAERO, &
864	ok_ade, ok_aie, flag_aerosol,flag_aerosol_strat) ! flags aerosols
865
866	! print *,'RADLWSW: apres RECMWF'
867	if(lldebug) then
868	CALL writefield_phy('zemtd_i',ZEMTD_i,klev+1)
869	CALL writefield_phy('zemtu_i',ZEMTU_i,klev+1)
870	CALL writefield_phy('ztrso_i',ZTRSO_i,klev+1)
871	CALL writefield_phy('zth_i',ZTH_i,klev+1)
872	CALL writefield_phy('zctrso',ZCTRSO,2)
873	CALL writefield_phy('zcemtr',ZCEMTR,2)
874	CALL writefield_phy('ztrsod',ZTRSOD,1)
875	CALL writefield_phy('zlwfc',ZLWFC,2)
876	CALL writefield_phy('zlwft_i',ZLWFT_i,klev+1)
877	CALL writefield_phy('zswfc',ZSWFC,2)
878	CALL writefield_phy('zswft_i',ZSWFT_i,klev+1)
879	CALL writefield_phy('psfswdir',PSFSWDIR,6)
880	CALL writefield_phy('psfswdif',PSFSWDIF,6)
881	CALL writefield_phy('pfsdnn',PFSDNN,1)
882	CALL writefield_phy('pfsdnv',PFSDNV,1)
883	CALL writefield_phy('ppiza_dst',PPIZA_TOT,klev)
884	CALL writefield_phy('pcga_dst',PCGA_TOT,klev)
885	CALL writefield_phy('ptaurel_dst',PTAU_TOT,klev)
886	CALL writefield_phy('zflux_i',ZFLUX_i,klev+1)
887	CALL writefield_phy('zfluc_i',ZFLUC_i,klev+1)
888	CALL writefield_phy('zfsdwn_i',ZFSDWN_i,klev+1)
889	CALL writefield_phy('zfsup_i',ZFSUP_i,klev+1)
890	CALL writefield_phy('zfcdwn_i',ZFCDWN_i,klev+1)
891	CALL writefield_phy('zfcup_i',ZFCUP_i,klev+1)
892	endif
893	! --------- output RECMWFL
894	! ZEMTD (KPROMA,KLEV+1) ; TOTAL DOWNWARD LONGWAVE EMISSIVITY
895	! ZEMTU (KPROMA,KLEV+1) ; TOTAL UPWARD LONGWAVE EMISSIVITY
896	! ZTRSO (KPROMA,KLEV+1) ; TOTAL SHORTWAVE TRANSMISSIVITY
897	! ZTH (KPROMA,KLEV+1) ; HALF LEVEL TEMPERATURE
898	! ZCTRSO (KPROMA,2) ; CLEAR-SKY SHORTWAVE TRANSMISSIVITY
899	! ZCEMTR (KPROMA,2) ; CLEAR-SKY NET LONGWAVE EMISSIVITY
900	! ZTRSOD (KPROMA) ; TOTAL-SKY SURFACE SW TRANSMISSITY
901	! ZLWFC (KPROMA,2) ; CLEAR-SKY LONGWAVE FLUXES
902	! ZLWFT (KPROMA,KLEV+1) ; TOTAL-SKY LONGWAVE FLUXES
903	! ZSWFC (KPROMA,2) ; CLEAR-SKY SHORTWAVE FLUXES
904	! ZSWFT (KPROMA,KLEV+1) ; TOTAL-SKY SHORTWAVE FLUXES
905	! PPIZA_TOT (KPROMA,KLEV,NSW); Single scattering albedo of total aerosols
906	! PCGA_TOT (KPROMA,KLEV,NSW); Assymetry factor for total aerosols
907	! PTAU_TOT (KPROMA,KLEV,NSW); Optical depth of total aerosols
908	! PPIZA_NAT (KPROMA,KLEV,NSW); Single scattering albedo of natural aerosols
909	! PCGA_NAT (KPROMA,KLEV,NSW); Assymetry factor for natural aerosols
910	! PTAU_NAT (KPROMA,KLEV,NSW); Optical depth of natiral aerosols
911	! PTAU_LW_TOT (KPROMA,KLEV,NLW); LW Optical depth of total aerosols
912	! PTAU_LW_NAT (KPROMA,KLEV,NLW); LW Optical depth of natural aerosols
913	! PSFSWDIR (KPROMA,NSW) ;
914	! PSFSWDIF (KPROMA,NSW) ;
915	! PFSDNN (KPROMA) ;
916	! PFSDNV (KPROMA) ;
917	! ---------
918	! ---------
919	! On retablit l'ordre des niveaux lmd pour les tableaux de sortie
920	! D autre part, on multiplie les resultats SW par fract pour etre coherent
921	! avec l ancien rayonnement AR4. Si nuit, fract=0 donc pas de
922	! rayonnement SW. (MPL 260609)
923	DO k=0,klev
924	DO i=1,klon
925	ZEMTD(i,k+1) = ZEMTD_i(i,k+1)
926	ZEMTU(i,k+1) = ZEMTU_i(i,k+1)
927	ZTRSO(i,k+1) = ZTRSO_i(i,k+1)
928	ZTH(i,k+1) = ZTH_i(i,k+1)
929	! ZLWFT(i,k+1) = ZLWFT_i(i,klev+1-k)
930	! ZSWFT(i,k+1) = ZSWFT_i(i,klev+1-k)
931	ZFLUP(i,k+1) = ZFLUX_i(i,1,k+1)
932	ZFLDN(i,k+1) = ZFLUX_i(i,2,k+1)
933	ZFLUP0(i,k+1) = ZFLUC_i(i,1,k+1)
934	ZFLDN0(i,k+1) = ZFLUC_i(i,2,k+1)
935	ZFSDN(i,k+1) = ZFSDWN_i(i,k+1)*fract(i)
936	ZFSDN0(i,k+1) = ZFCDWN_i(i,k+1)*fract(i)
937	ZFSUP (i,k+1) = ZFSUP_i(i,k+1)*fract(i)
938	ZFSUP0(i,k+1) = ZFCUP_i(i,k+1)*fract(i)
939	! Nouveau calcul car visiblement ZSWFT et ZSWFC sont nuls dans RRTM cy32
940	! en sortie de radlsw.F90 - MPL 7.01.09
941	ZSWFT(i,k+1) = (ZFSDWN_i(i,k+1)-ZFSUP_i(i,k+1))*fract(i)
942	ZSWFT0_i(i,k+1) = (ZFCDWN_i(i,k+1)-ZFCUP_i(i,k+1))*fract(i)
943	! WRITE(*,'("FSDN FSUP FCDN FCUP: ",4E12.5)') ZFSDWN_i(i,k+1),&
944	! ZFSUP_i(i,k+1),ZFCDWN_i(i,k+1),ZFCUP_i(i,k+1)
945	ZLWFT(i,k+1) =-ZFLUX_i(i,2,k+1)-ZFLUX_i(i,1,k+1)
946	ZLWFT0_i(i,k+1)=-ZFLUC_i(i,2,k+1)-ZFLUC_i(i,1,k+1)
947	! print *,'FLUX2 FLUX1 FLUC2 FLUC1',ZFLUX_i(i,2,k+1),&
948	! & ZFLUX_i(i,1,k+1),ZFLUC_i(i,2,k+1),ZFLUC_i(i,1,k+1)
949	ENDDO
950	ENDDO
951
952	!--ajout OB
953	ZTOPSWADAERO(:) =ZTOPSWADAERO(:) *fract(:)
954	ZSOLSWADAERO(:) =ZSOLSWADAERO(:) *fract(:)
955	ZTOPSWAD0AERO(:)=ZTOPSWAD0AERO(:)*fract(:)
956	ZSOLSWAD0AERO(:)=ZSOLSWAD0AERO(:)*fract(:)
957	ZTOPSWAIAERO(:) =ZTOPSWAIAERO(:) *fract(:)
958	ZSOLSWAIAERO(:) =ZSOLSWAIAERO(:) *fract(:)
959	ZTOPSWCF_AERO(:,1)=ZTOPSWCF_AERO(:,1)*fract(:)
960	ZTOPSWCF_AERO(:,2)=ZTOPSWCF_AERO(:,2)*fract(:)
961	ZTOPSWCF_AERO(:,3)=ZTOPSWCF_AERO(:,3)*fract(:)
962	ZSOLSWCF_AERO(:,1)=ZSOLSWCF_AERO(:,1)*fract(:)
963	ZSOLSWCF_AERO(:,2)=ZSOLSWCF_AERO(:,2)*fract(:)
964	ZSOLSWCF_AERO(:,3)=ZSOLSWCF_AERO(:,3)*fract(:)
965
966	! print*,'SW_RRTM ZFSDN0 1 , klev:',ZFSDN0(1:klon,1),ZFSDN0(1:klon,klev)
967	! print*,'SW_RRTM ZFSUP0 1 , klev:',ZFSUP0(1:klon,1),ZFSUP0(1:klon,klev)
968	! print*,'SW_RRTM ZFSDN 1 , klev:',ZFSDN(1:klon,1),ZFSDN(1:klon,klev)
969	! print*,'SW_RRTM ZFSUP 1 , klev:',ZFSUP(1:klon,1),ZFSUP(1:klon,klev)
970	! print*,'OK1'
971	! ---------
972	! ---------
973	! On renseigne les champs LMDz, pour avoir la meme chose qu'en sortie de
974	! LW_LMDAR4 et SW_LMDAR4
975	DO i = 1, kdlon
976	zsolsw(i) = ZSWFT(i,1)
977	zsolsw0(i) = ZSWFT0_i(i,1)
978	! zsolsw0(i) = ZFSDN0(i,1) -ZFSUP0(i,1)
979	ztopsw(i) = ZSWFT(i,klev+1)
980	ztopsw0(i) = ZSWFT0_i(i,klev+1)
981	! ztopsw0(i) = ZFSDN0(i,klev+1)-ZFSUP0(i,klev+1)
982	!
983	! zsollw(i) = ZFLDN(i,1) -ZFLUP(i,1)
984	! zsollw0(i) = ZFLDN0(i,1) -ZFLUP0(i,1)
985	! ztoplw(i) = ZFLDN(i,klev+1) -ZFLUP(i,klev+1)
986	! ztoplw0(i) = ZFLDN0(i,klev+1)-ZFLUP0(i,klev+1)
987	zsollw(i) = ZLWFT(i,1)
988	zsollw0(i) = ZLWFT0_i(i,1)
989	ztoplw(i) = ZLWFT(i,klev+1)*(-1)
990	ztoplw0(i) = ZLWFT0_i(i,klev+1)*(-1)
991	!
992	IF (fract(i) == 0.) THEN
993	!!!!! A REVOIR MPL (20090630) ca n a pas de sens quand fract=0
994	! pas plus que dans le sw_AR4
995	zalbpla(i) = 1.0e+39
996	ELSE
997	zalbpla(i) = ZFSUP(i,klev+1)/ZFSDN(i,klev+1)
998	ENDIF
999	!!! 5 juin 2015
1000	!!! Correction MP bug RRTM
1001	zsollwdown(i)= -1.*ZFLDN(i,1)
1002	ENDDO
1003	! print*,'OK2'
1004
1005	! extrait de SW_AR4
1006	! DO k = 1, KFLEV
1007	! kpl1 = k+1
1008	! DO i = 1, KDLON
1009	! PHEAT(i,k) = -(ZFSUP(i,kpl1)-ZFSUP(i,k)) -(ZFSDN(i,k)-ZFSDN(i,kpl1))
1010	! PHEAT(i,k) = PHEAT(i,k) * RDAY*RG/RCPD / PDP(i,k)
1011	! ZLWFT(klon,k),ZSWFT
1012
1013	do k=1,kflev
1014	do i=1,kdlon
1015	zheat(i,k)=(ZSWFT(i,k+1)-ZSWFT(i,k))RDAYRG/RCPD/PDP(i,k)
1016	zheat0(i,k)=(ZSWFT0_i(i,k+1)-ZSWFT0_i(i,k))RDAYRG/RCPD/PDP(i,k)
1017	zcool(i,k)=(ZLWFT(i,k)-ZLWFT(i,k+1))RDAYRG/RCPD/PDP(i,k)
1018	zcool0(i,k)=(ZLWFT0_i(i,k)-ZLWFT0_i(i,k+1))RDAYRG/RCPD/PDP(i,k)
1019	! print *,'heat cool heat0 cool0 ',zheat(i,k),zcool(i,k),zheat0(i,k),zcool0(i,k)
1020	! ZFLUCUP_i(i,k)=ZFLUC_i(i,1,k)
1021	! ZFLUCDWN_i(i,k)=ZFLUC_i(i,2,k)
1022	enddo
1023	enddo
1024	#else
1025	abort_message="You should compile with -rrtm if running with iflag_rrtm=1"
1026	call abort_physic(modname, abort_message, 1)
1027	#endif
1028	ENDIF ! iflag_rrtm
1029	!======================================================================
1030
1031	DO i = 1, kdlon
1032	topsw(iof+i) = ztopsw(i)
1033	toplw(iof+i) = ztoplw(i)
1034	solsw(iof+i) = zsolsw(i)
1035	sollw(iof+i) = zsollw(i)
1036	sollwdown(iof+i) = zsollwdown(i)
1037	DO k = 1, kflev+1
1038	lwdn0 ( iof+i,k) = ZFLDN0 ( i,k)
1039	lwdn ( iof+i,k) = ZFLDN ( i,k)
1040	lwup0 ( iof+i,k) = ZFLUP0 ( i,k)
1041	lwup ( iof+i,k) = ZFLUP ( i,k)
1042	ENDDO
1043	topsw0(iof+i) = ztopsw0(i)
1044	toplw0(iof+i) = ztoplw0(i)
1045	solsw0(iof+i) = zsolsw0(i)
1046	sollw0(iof+i) = zsollw0(i)
1047	albpla(iof+i) = zalbpla(i)
1048
1049	DO k = 1, kflev+1
1050	swdn0 ( iof+i,k) = ZFSDN0 ( i,k)
1051	swdn ( iof+i,k) = ZFSDN ( i,k)
1052	swup0 ( iof+i,k) = ZFSUP0 ( i,k)
1053	swup ( iof+i,k) = ZFSUP ( i,k)
1054	ENDDO
1055	ENDDO
1056	!-transform the aerosol forcings, if they have
1057	! to be calculated
1058	IF (ok_ade) THEN
1059	DO i = 1, kdlon
1060	topswad_aero(iof+i) = ztopswadaero(i)
1061	topswad0_aero(iof+i) = ztopswad0aero(i)
1062	solswad_aero(iof+i) = zsolswadaero(i)
1063	solswad0_aero(iof+i) = zsolswad0aero(i)
1064	! MS the following lines seem to be wrong, why is iof on right hand side???
1065	! topsw_aero(iof+i,:) = ztopsw_aero(iof+i,:)
1066	! topsw0_aero(iof+i,:) = ztopsw0_aero(iof+i,:)
1067	! solsw_aero(iof+i,:) = zsolsw_aero(iof+i,:)
1068	! solsw0_aero(iof+i,:) = zsolsw0_aero(iof+i,:)
1069	topsw_aero(iof+i,:) = ztopsw_aero(i,:)
1070	topsw0_aero(iof+i,:) = ztopsw0_aero(i,:)
1071	solsw_aero(iof+i,:) = zsolsw_aero(i,:)
1072	solsw0_aero(iof+i,:) = zsolsw0_aero(i,:)
1073	topswcf_aero(iof+i,:) = ztopswcf_aero(i,:)
1074	solswcf_aero(iof+i,:) = zsolswcf_aero(i,:)
1075	!-LW
1076	toplwad_aero(iof+i) = ztoplwadaero(i)
1077	toplwad0_aero(iof+i) = ztoplwad0aero(i)
1078	sollwad_aero(iof+i) = zsollwadaero(i)
1079	sollwad0_aero(iof+i) = zsollwad0aero(i)
1080	ENDDO
1081	ELSE
1082	DO i = 1, kdlon
1083	topswad_aero(iof+i) = 0.0
1084	solswad_aero(iof+i) = 0.0
1085	topswad0_aero(iof+i) = 0.0
1086	solswad0_aero(iof+i) = 0.0
1087	topsw_aero(iof+i,:) = 0.
1088	topsw0_aero(iof+i,:) =0.
1089	solsw_aero(iof+i,:) = 0.
1090	solsw0_aero(iof+i,:) = 0.
1091	!-LW
1092	toplwad_aero(iof+i) = 0.0
1093	sollwad_aero(iof+i) = 0.0
1094	toplwad0_aero(iof+i) = 0.0
1095	sollwad0_aero(iof+i) = 0.0
1096	ENDDO
1097	ENDIF
1098	IF (ok_aie) THEN
1099	DO i = 1, kdlon
1100	topswai_aero(iof+i) = ztopswaiaero(i)
1101	solswai_aero(iof+i) = zsolswaiaero(i)
1102	!-LW
1103	toplwai_aero(iof+i) = ztoplwaiaero(i)
1104	sollwai_aero(iof+i) = zsollwaiaero(i)
1105	ENDDO
1106	ELSE
1107	DO i = 1, kdlon
1108	topswai_aero(iof+i) = 0.0
1109	solswai_aero(iof+i) = 0.0
1110	!-LW
1111	toplwai_aero(iof+i) = 0.0
1112	sollwai_aero(iof+i) = 0.0
1113	ENDDO
1114	ENDIF
1115	DO k = 1, kflev
1116	DO i = 1, kdlon
1117	! scale factor to take into account the difference between
1118	! dry air and watter vapour scpecifi! heat capacity
1119	zznormcp=1.0+RVTMP2*PWV(i,k)
1120	heat(iof+i,k) = zheat(i,k)/zznormcp
1121	cool(iof+i,k) = zcool(i,k)/zznormcp
1122	heat0(iof+i,k) = zheat0(i,k)/zznormcp
1123	cool0(iof+i,k) = zcool0(i,k)/zznormcp
1124	ENDDO
1125	ENDDO
1126
1127	ENDDO ! j = 1, nb_gr
1128
1129	END SUBROUTINE radlwsw
1130
1131	end module radlwsw_m

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

Download in other formats:

Original Format