Context Navigation

radlwsw_m.F90

Last change on this file was 2003, checked in by Laurent Fairhead, 11 years ago

Nouvelle version qui inclut les effets des aérosols et propose les mêmes diagnostics des effets
directs et indirects que l'ancienne version du rayonnement.
OB

New RRTM version that includes the effects of aerosols and outputs the same direct and indirect effects
diagnostics as the old version
OB

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: 41.2 KB

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

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

Download in other formats:

Original Format