Context Navigation

radlwsw_m.F90 @ 3817

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

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

Download in other formats:

Original Format