Context Navigation

radlwsw.F @ 1134

Last change on this file since 1134 was 998, checked in by Laurent Fairhead, 16 years ago
Modifications necessaires a la preparation au passage au nouveau rayonnement RRTM MPL LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 14.7 KB

Line
1	!
2	! $Header$
3	!
4	SUBROUTINE radlwsw(dist, rmu0, fract,
5	. paprs, pplay,tsol,alb1, alb2, t,q,wo,
6	. cldfra, cldemi, cldtaupd,
7	. heat,heat0,cool,cool0,radsol,albpla,
8	. topsw,toplw,solsw,sollw,
9	. sollwdown,
10	. topsw0,toplw0,solsw0,sollw0,
11	. lwdn0, lwdn, lwup0, lwup,
12	. swdn0, swdn, swup0, swup,
13	. ok_ade, ok_aie,
14	. tau_ae, piz_ae, cg_ae,
15	. topswad, solswad,
16	. cldtaupi, topswai, solswai,qsat,flwc,fiwc)
17	c
18	USE dimphy
19	IMPLICIT none
20	c======================================================================
21	c Auteur(s): Z.X. Li (LMD/CNRS) date: 19960719
22	c Objet: interface entre le modele et les rayonnements
23	c Arguments:
24	c dist-----input-R- distance astronomique terre-soleil
25	c rmu0-----input-R- cosinus de l'angle zenithal
26	c fract----input-R- duree d'ensoleillement normalisee
27	c co2_ppm--input-R- concentration du gaz carbonique (en ppm)
28	c solaire--input-R- constante solaire (W/m**2)
29	c paprs----input-R- pression a inter-couche (Pa)
30	c pplay----input-R- pression au milieu de couche (Pa)
31	c tsol-----input-R- temperature du sol (en K)
32	c alb1-----input-R- albedo du sol(entre 0 et 1) dans l'interval visible
33	c alb2-----input-R- albedo du sol(entre 0 et 1) dans l'interval proche infra-rouge
34	c t--------input-R- temperature (K)
35	c q--------input-R- vapeur d'eau (en kg/kg)
36	c wo-------input-R- contenu en ozone (en kg/kg) correction MPL 100505
37	c cldfra---input-R- fraction nuageuse (entre 0 et 1)
38	c cldtaupd---input-R- epaisseur optique des nuages dans le visible (present-day value)
39	c cldemi---input-R- emissivite des nuages dans l'IR (entre 0 et 1)
40	c ok_ade---input-L- apply the Aerosol Direct Effect or not?
41	c ok_aie---input-L- apply the Aerosol Indirect Effect or not?
42	c tau_ae, piz_ae, cg_ae-input-R- aerosol optical properties (calculated in aeropt.F)
43	c cldtaupi-input-R- epaisseur optique des nuages dans le visible
44	c calculated for pre-industrial (pi) aerosol concentrations, i.e. with smaller
45	c droplet concentration, thus larger droplets, thus generally cdltaupi cldtaupd
46	c it is needed for the diagnostics of the aerosol indirect radiative forcing
47	c
48	c heat-----output-R- echauffement atmospherique (visible) (K/jour)
49	c cool-----output-R- refroidissement dans l'IR (K/jour)
50	c radsol---output-R- bilan radiatif net au sol (W/m**2) (+ vers le bas)
51	c albpla---output-R- albedo planetaire (entre 0 et 1)
52	c topsw----output-R- flux solaire net au sommet de l'atm.
53	c toplw----output-R- ray. IR montant au sommet de l'atmosphere
54	c solsw----output-R- flux solaire net a la surface
55	c sollw----output-R- ray. IR montant a la surface
56	c solswad---output-R- ray. solaire net absorbe a la surface (aerosol dir)
57	c topswad---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol dir)
58	c solswai---output-R- ray. solaire net absorbe a la surface (aerosol ind)
59	c topswai---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol ind)
60	c
61	c ATTENTION: swai and swad have to be interpreted in the following manner:
62	c ---------
63	c ok_ade=F & ok_aie=F -both are zero
64	c ok_ade=T & ok_aie=F -aerosol direct forcing is F_{AD} = topsw-topswad
65	c indirect is zero
66	c ok_ade=F & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai
67	c direct is zero
68	c ok_ade=T & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai
69	c aerosol direct forcing is F_{AD} = topswai-topswad
70	c
71
72	c======================================================================
73	cym#include "dimensions.h"
74	cym#include "dimphy.h"
75	cym#include "raddim.h"
76	#include "YOETHF.h"
77	c
78	real rmu0(klon), fract(klon), dist
79	cIM real co2_ppm
80	cIM real solaire
81	#include "clesphys.h"
82	c
83	real paprs(klon,klev+1), pplay(klon,klev)
84	real alb1(klon), alb2(klon), tsol(klon)
85	real t(klon,klev), q(klon,klev), wo(klon,klev)
86	real cldfra(klon,klev), cldemi(klon,klev), cldtaupd(klon,klev)
87	real heat(klon,klev), cool(klon,klev)
88	real heat0(klon,klev), cool0(klon,klev)
89	real radsol(klon), topsw(klon), toplw(klon)
90	real solsw(klon), sollw(klon), albpla(klon)
91	real topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)
92	real sollwdown(klon)
93	cIM output 3D
94	REAL*8 ZFSUP(KDLON,KFLEV+1)
95	REAL*8 ZFSDN(KDLON,KFLEV+1)
96	REAL*8 ZFSUP0(KDLON,KFLEV+1)
97	REAL*8 ZFSDN0(KDLON,KFLEV+1)
98	c
99	REAL*8 ZFLUP(KDLON,KFLEV+1)
100	REAL*8 ZFLDN(KDLON,KFLEV+1)
101	REAL*8 ZFLUP0(KDLON,KFLEV+1)
102	REAL*8 ZFLDN0(KDLON,KFLEV+1)
103	c
104	REAL*8 zx_alpha1, zx_alpha2
105	c
106	#include "YOMCST.h"
107	c
108	INTEGER k, kk, i, j, iof, nb_gr
109	EXTERNAL LW_LMDAR4,SW_LMDAR4
110	c
111	cIM ctes ds clesphys.h REAL*8 RCO2, RCH4, RN2O, RCFC11, RCFC12
112	REAL*8 PSCT
113	c
114	REAL*8 PALBD(kdlon,2), PALBP(kdlon,2)
115	REAL*8 PEMIS(kdlon), PDT0(kdlon), PVIEW(kdlon)
116	REAL*8 PPSOL(kdlon), PDP(kdlon,klev)
117	REAL*8 PTL(kdlon,kflev+1), PPMB(kdlon,kflev+1)
118	REAL*8 PTAVE(kdlon,kflev)
119	REAL*8 PWV(kdlon,kflev), PQS(kdlon,kflev), POZON(kdlon,kflev)
120	REAL*8 PAER(kdlon,kflev,5)
121	REAL*8 PCLDLD(kdlon,kflev)
122	REAL*8 PCLDLU(kdlon,kflev)
123	REAL*8 PCLDSW(kdlon,kflev)
124	REAL*8 PTAU(kdlon,2,kflev)
125	REAL*8 POMEGA(kdlon,2,kflev)
126	REAL*8 PCG(kdlon,2,kflev)
127	c
128	REAL*8 zfract(kdlon), zrmu0(kdlon), zdist
129	c
130	REAL*8 zheat(kdlon,kflev), zcool(kdlon,kflev)
131	REAL*8 zheat0(kdlon,kflev), zcool0(kdlon,kflev)
132	REAL*8 ztopsw(kdlon), ztoplw(kdlon)
133	REAL*8 zsolsw(kdlon), zsollw(kdlon), zalbpla(kdlon)
134	cIM
135	REAL*8 zsollwdown(kdlon)
136	c
137	REAL*8 ztopsw0(kdlon), ztoplw0(kdlon)
138	REAL*8 zsolsw0(kdlon), zsollw0(kdlon)
139	REAL*8 zznormcp
140	cIM output 3D : SWup, SWdn, LWup, LWdn
141	REAL swdn(klon,kflev+1),swdn0(klon,kflev+1)
142	REAL swup(klon,kflev+1),swup0(klon,kflev+1)
143	REAL lwdn(klon,kflev+1),lwdn0(klon,kflev+1)
144	REAL lwup(klon,kflev+1),lwup0(klon,kflev+1)
145	REAL qsat(klon,klev),flwc(klon,klev),fiwc(klon,klev)
146	c-OB
147	cjq the following quantities are needed for the aerosol radiative forcings
148
149	real topswad(klon), solswad(klon) ! output: aerosol direct forcing at TOA and surface
150	real topswai(klon), solswai(klon) ! output: aerosol indirect forcing atTOA and surface
151	real tau_ae(klon,klev,2), piz_ae(klon,klev,2), cg_ae(klon,klev,2) ! aerosol optical properties (see aeropt.F)
152	real cldtaupi(klon,klev) ! cloud optical thickness for pre-industrial aerosol concentrations
153	! (i.e., with a smaller droplet concentrationand thus larger droplet radii)
154	logical ok_ade, ok_aie ! switches whether to use aerosol direct (indirect) effects or not
155	real*8 tauae(kdlon,kflev,2) ! aer opt properties
156	real*8 pizae(kdlon,kflev,2)
157	real*8 cgae(kdlon,kflev,2)
158	REAL*8 PTAUA(kdlon,2,kflev) ! present-day value of cloud opt thickness (PTAU is pre-industrial value), local use
159	REAL*8 POMEGAA(kdlon,2,kflev) ! dito for single scatt albedo
160	REAL*8 ztopswad(kdlon), zsolswad(kdlon) ! Aerosol direct forcing at TOAand surface
161	REAL*8 ztopswai(kdlon), zsolswai(kdlon) ! dito, indirect
162	cjq-end
163	!rv
164	tauae(:,:,:)=0.
165	pizae(:,:,:)=0.
166	cgae(:,:,:)=0.
167	!rv
168
169	c
170	c-------------------------------------------
171	nb_gr = klon / kdlon
172	IF (nb_gr*kdlon .NE. klon) THEN
173	PRINT*, "kdlon mauvais:", klon, kdlon, nb_gr
174	CALL abort
175	ENDIF
176	IF (kflev .NE. klev) THEN
177	PRINT*, "kflev differe de klev, kflev, klev"
178	CALL abort
179	ENDIF
180	c-------------------------------------------
181	DO k = 1, klev
182	DO i = 1, klon
183	heat(i,k)=0.
184	cool(i,k)=0.
185	heat0(i,k)=0.
186	cool0(i,k)=0.
187	ENDDO
188	ENDDO
189	c
190	zdist = dist
191	c
192	cIM anciennes valeurs
193	c RCO2 = co2_ppm * 1.0e-06 * 44.011/28.97
194	c
195	cIM : on met RCO2, RCH4, RN2O, RCFC11 et RCFC12 dans clesphys.h /lecture ds conf_phys.F90
196	c RCH4 = 1.65E-06* 16.043/28.97
197	c RN2O = 306.E-09* 44.013/28.97
198	c RCFC11 = 280.E-12* 137.3686/28.97
199	c RCFC12 = 484.E-12* 120.9140/28.97
200	cIM anciennes valeurs
201	c RCH4 = 1.72E-06* 16.043/28.97
202	c RN2O = 310.E-09* 44.013/28.97
203	c
204	c PRINT*,'IMradlwsw : solaire, co2= ', solaire, co2_ppm
205	PSCT = solaire/zdist/zdist
206	c
207	DO 99999 j = 1, nb_gr
208	iof = kdlon*(j-1)
209	c
210	DO i = 1, kdlon
211	zfract(i) = fract(iof+i)
212	zrmu0(i) = rmu0(iof+i)
213	PALBD(i,1) = alb1(iof+i)
214	! PALBD(i,2) = alb1(iof+i)
215	PALBD(i,2) = alb2(iof+i)
216	PALBP(i,1) = alb1(iof+i)
217	! PALBP(i,2) = alb1(iof+i)
218	PALBP(i,2) = alb2(iof+i)
219	cIM cf. JLD pour etre en accord avec ORCHIDEE il faut mettre PEMIS(i) = 0.96
220	PEMIS(i) = 1.0
221	PVIEW(i) = 1.66
222	PPSOL(i) = paprs(iof+i,1)
223	zx_alpha1 = (paprs(iof+i,1)-pplay(iof+i,2))
224	. / (pplay(iof+i,1)-pplay(iof+i,2))
225	zx_alpha2 = 1.0 - zx_alpha1
226	PTL(i,1) = t(iof+i,1) * zx_alpha1 + t(iof+i,2) * zx_alpha2
227	PTL(i,klev+1) = t(iof+i,klev)
228	PDT0(i) = tsol(iof+i) - PTL(i,1)
229	ENDDO
230	DO k = 2, kflev
231	DO i = 1, kdlon
232	PTL(i,k) = (t(iof+i,k)+t(iof+i,k-1))*0.5
233	ENDDO
234	ENDDO
235	DO k = 1, kflev
236	DO i = 1, kdlon
237	PDP(i,k) = paprs(iof+i,k)-paprs(iof+i,k+1)
238	PTAVE(i,k) = t(iof+i,k)
239	PWV(i,k) = MAX (q(iof+i,k), 1.0e-12)
240	PQS(i,k) = PWV(i,k)
241	c wo: cm.atm (epaisseur en cm dans la situation standard)
242	c POZON: kg/kg
243	POZON(i,k) = MAX(wo(iof+i,k),1.0e-12)*RG/46.6968
244	. /(paprs(iof+i,k)-paprs(iof+i,k+1))
245	. *(paprs(iof+i,1)/101325.0)
246	PCLDLD(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)
247	PCLDLU(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)
248	PCLDSW(i,k) = cldfra(iof+i,k)
249	PTAU(i,1,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! 1e-12 serait instable
250	PTAU(i,2,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! pour 32-bit machines
251	POMEGA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAU(i,1,k))
252	POMEGA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAU(i,2,k))
253	PCG(i,1,k) = 0.865
254	PCG(i,2,k) = 0.910
255	c-OB
256	cjq Introduced for aerosol indirect forcings.
257	cjq The following values use the cloud optical thickness calculated from
258	cjq present-day aerosol concentrations whereas the quantities without the
259	cjq "A" at the end are for pre-industial (natural-only) aerosol concentrations
260	cjq
261	PTAUA(i,1,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! 1e-12 serait instable
262	PTAUA(i,2,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! pour 32-bit machines
263	POMEGAA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAUA(i,1,k))
264	POMEGAA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAUA(i,2,k))
265	cjq-end
266	ENDDO
267	ENDDO
268	c
269	DO k = 1, kflev+1
270	DO i = 1, kdlon
271	PPMB(i,k) = paprs(iof+i,k)/100.0
272	ENDDO
273	ENDDO
274	c
275	DO kk = 1, 5
276	DO k = 1, kflev
277	DO i = 1, kdlon
278	PAER(i,k,kk) = 1.0E-15
279	ENDDO
280	ENDDO
281	ENDDO
282	c-OB
283	DO k = 1, kflev
284	DO i = 1, kdlon
285	tauae(i,k,1)=tau_ae(iof+i,k,1)
286	pizae(i,k,1)=piz_ae(iof+i,k,1)
287	cgae(i,k,1) =cg_ae(iof+i,k,1)
288	tauae(i,k,2)=tau_ae(iof+i,k,2)
289	pizae(i,k,2)=piz_ae(iof+i,k,2)
290	cgae(i,k,2) =cg_ae(iof+i,k,2)
291	ENDDO
292	ENDDO
293	c
294	c===== si iflag_rrtm=0 ================================================
295	cIM ctes ds clesphys.h CALL LW(RCO2,RCH4,RN2O,RCFC11,RCFC12,
296	cIM ctes ds clesphys.h CALL SW(PSCT, RCO2, zrmu0, zfract,
297	c
298	if (iflag_rrtm.eq.0) then
299	CALL LW_LMDAR4(
300	. PPMB, PDP,
301	. PPSOL,PDT0,PEMIS,
302	. PTL, PTAVE, PWV, POZON, PAER,
303	. PCLDLD,PCLDLU,
304	. PVIEW,
305	. zcool, zcool0,
306	. ztoplw,zsollw,ztoplw0,zsollw0,
307	. zsollwdown,
308	. ZFLUP, ZFLDN, ZFLUP0,ZFLDN0)
309	CALL SW_LMDAR4(PSCT, zrmu0, zfract,
310	S PPMB, PDP,
311	S PPSOL, PALBD, PALBP,
312	S PTAVE, PWV, PQS, POZON, PAER,
313	S PCLDSW, PTAU, POMEGA, PCG,
314	S zheat, zheat0,
315	S zalbpla,ztopsw,zsolsw,ztopsw0,zsolsw0,
316	S ZFSUP,ZFSDN,ZFSUP0,ZFSDN0,
317	S tauae, pizae, cgae, ! aerosol optical properties
318	s PTAUA, POMEGAA,
319	s ztopswad,zsolswad,ztopswai,zsolswai, ! diagnosed aerosol forcing
320	J ok_ade, ok_aie) ! apply aerosol effects or not?
321	else
322	c===== si iflag_rrtm=1, on passe dans SW via RECMWFL ===============
323	PRINT*, "Cette option ne fonctionne pas encore !!!"
324	CALL abort
325	endif ! if(iflag_rrtm=0)
326
327	c======================================================================
328	DO i = 1, kdlon
329	radsol(iof+i) = zsolsw(i) + zsollw(i)
330	topsw(iof+i) = ztopsw(i)
331	toplw(iof+i) = ztoplw(i)
332	solsw(iof+i) = zsolsw(i)
333	sollw(iof+i) = zsollw(i)
334	sollwdown(iof+i) = zsollwdown(i)
335	cIM
336	DO k = 1, kflev+1
337	lwdn0 ( iof+i,k) = ZFLDN0 ( i,k)
338	lwdn ( iof+i,k) = ZFLDN ( i,k)
339	lwup0 ( iof+i,k) = ZFLUP0 ( i,k)
340	lwup ( iof+i,k) = ZFLUP ( i,k)
341	ENDDO
342	c
343	topsw0(iof+i) = ztopsw0(i)
344	toplw0(iof+i) = ztoplw0(i)
345	solsw0(iof+i) = zsolsw0(i)
346	sollw0(iof+i) = zsollw0(i)
347	albpla(iof+i) = zalbpla(i)
348	cIM
349	DO k = 1, kflev+1
350	swdn0 ( iof+i,k) = ZFSDN0 ( i,k)
351	swdn ( iof+i,k) = ZFSDN ( i,k)
352	swup0 ( iof+i,k) = ZFSUP0 ( i,k)
353	swup ( iof+i,k) = ZFSUP ( i,k)
354	ENDDO !k=1, kflev+1
355	ENDDO
356	cjq-transform the aerosol forcings, if they have
357	cjq to be calculated
358	IF (ok_ade) THEN
359	DO i = 1, kdlon
360	topswad(iof+i) = ztopswad(i)
361	solswad(iof+i) = zsolswad(i)
362	ENDDO
363	ELSE
364	DO i = 1, kdlon
365	topswad(iof+i) = 0.0
366	solswad(iof+i) = 0.0
367	ENDDO
368	ENDIF
369	IF (ok_aie) THEN
370	DO i = 1, kdlon
371	topswai(iof+i) = ztopswai(i)
372	solswai(iof+i) = zsolswai(i)
373	ENDDO
374	ELSE
375	DO i = 1, kdlon
376	topswai(iof+i) = 0.0
377	solswai(iof+i) = 0.0
378	ENDDO
379	ENDIF
380	cjq-end
381	DO k = 1, kflev
382	c DO i = 1, kdlon
383	c heat(iof+i,k) = zheat(i,k)
384	c cool(iof+i,k) = zcool(i,k)
385	c heat0(iof+i,k) = zheat0(i,k)
386	c cool0(iof+i,k) = zcool0(i,k)
387	c ENDDO
388	DO i = 1, kdlon
389	C scale factor to take into account the difference between
390	C dry air and watter vapour scpecific heat capacity
391	zznormcp=1.0+RVTMP2*PWV(i,k)
392	heat(iof+i,k) = zheat(i,k)/zznormcp
393	cool(iof+i,k) = zcool(i,k)/zznormcp
394	heat0(iof+i,k) = zheat0(i,k)/zznormcp
395	cool0(iof+i,k) = zcool0(i,k)/zznormcp
396	ENDDO
397	ENDDO
398	c
399	99999 CONTINUE
400	RETURN
401	END

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

Download in other formats:

Original Format