Context Navigation

radlwsw_m.F90 @ 3814

Last change on this file since 3814 was 3814, checked in by ymipsl, 10 years ago

remove all dynamic dependency in LMDZ physics except for the include "dimensions.h"

File size: 45.6 KB

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

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

Download in other formats:

Original Format