Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

rforcing.F @ 3400

Last change on this file since 3400 was 1547, checked in by idelkadi, 13 years ago
Rajout des routines pour le calcul du forcage radiatif "offline" dans LMDZ ! Note technique disponible sur : http://www.lmd.jussieu.fr/~idelkadi/Welcome_diagnos.html
Property svn:executable set to ``*
File size: 50.7 KB

Line
1	C%%%%%%%%%%%%%%%%%% Abderrahmane 01 06 2005 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2
3	SUBROUTINE RFORCING ( m,rjour,read_oz,
4	e tfi,qfi,wo,wo_dl,wo_ref,wo_dl_ref,
5	e pplayfi,paprsfi,cldfrafi,cldliqfi,
6	e tsolfi,albsfi,rlatfi,
7	s t_seri,
8	s tps_sw_ref,tps_sw_ref0,tps_lw_ref,tps_lw_ref0,
9	s tps_sw_ini,tps_sw_ini0,tps_lw_ini,tps_lw_ini0,
10	s d_tps_sw_ini,d_tps_sw_ini0,d_tps_lw_ini,d_tps_lw_ini0,
11	s d_tps_sw_adj,d_tps_sw_adj0,d_tps_lw_adj,d_tps_lw_adj0,
12	s toa_sw_ref,toa_sw_ref0,toa_lw_ref,toa_lw_ref0,
13	s toa_sw_ini,toa_sw_ini0,toa_lw_ini,toa_lw_ini0,
14	s d_toa_sw_ini,d_toa_sw_ini0,d_toa_lw_ini,d_toa_lw_ini0,
15	s d_toa_sw_adj,d_toa_sw_adj0,d_toa_lw_adj,d_toa_lw_adj0,
16	s srf_sw_ref,srf_sw_ref0,srf_lw_ref,srf_lw_ref0,
17	s d_srf_sw_adj,d_srf_sw_adj0,d_srf_lw_adj,d_srf_lw_adj0,
18	s dHrad_dT,bilq_ref,bilq_ini,bilq_adj,
19	s heat_ref, heat0_ref, cool_ref, cool0_ref,
20	s d_heat_ini, d_heat0_ini, d_cool_ini, d_cool0_ini,
21	s d_heat_adj, d_heat0_adj, d_cool_adj, d_cool0_adj )
22	C
23	C output
24	C t_seri: adjusted temperature profile
25	!
26	! JLD: 2009-07-13
27	! Improvment of the convergence of the temperature adjustment.
28	! Variable dHrad_dT has been introduced.
29	!
30	! JLD: 2009-07-12
31	! add the flux and the forcing at the surface
32	C
33	use dimphy
34	use IOIPSL
35	! use ozonecm_m, only: ozonecm ! ozone of J.-F. Royer
36	use radlwsw_m, only: radlwsw
37
38	IMPLICIT none
39
40	#include "clesphys.h"
41	#include "nuage.h"
42	#include "thermcell.h"
43	#include "YOMCST.h"
44
45	LOGICAL debug
46	PARAMETER (debug=.false.)
47	REAL rjour,dtime
48	PARAMETER (dtime=4.*86400.0)
49	INTEGER i, k, j, itap, m, nbs
50	PARAMETER (nbs=30*2)
51	CJLD PARAMETER (nbs=2)
52
53	REAL co2_ppm_ref, CH4_ppb_ref, N2O_ppb_ref
54	SAVE co2_ppm_ref, CH4_ppb_ref, N2O_ppb_ref
55
56	REAL CFC11_ppt_ref, CFC12_ppt_ref
57	SAVE CFC11_ppt_ref, CFC12_ppt_ref
58
59	REAL rco2_ref
60	REAL rch4_ref
61	REAL RN2O_ref
62	Real rcfc11_ref
63	Real rcfc12_ref
64	REAL solaire_ref
65	SAVE rco2_ref,rch4_ref,RN2O_ref,
66	$ rcfc11_ref,rcfc12_ref,solaire_ref
67
68	REAL co2_ppm_modif,solaire_modif
69	REAL RCO2_modif,RCH4_modif,RN2O_modif,
70	$ RCFC11_modif,RCFC12_modif
71	REAL CH4_ppb_modif,N2O_ppb_modif,
72	$ CFC11_ppt_modif,CFC12_ppt_modif
73	SAVE RCO2_modif,RCH4_modif,RN2O_modif,
74	$ RCFC11_modif,RCFC12_modif,solaire_modif
75
76	c initialisation des constantes
77
78	logical ok_newmicro
79	data ok_newmicro/.true./
80	save ok_newmicro
81
82
83	REAL pct(klon), pctlwp(klon), pch(klon), pcl(klon), pcm(klon)
84	REAL xflwp(klon), xfiwp(klon)
85	REAL xflwc(klon,klev), xfiwc(klon,klev)
86	REAL sulfate_ref(klon, klev) ! sulfate aerosol mass concentration [ug m-3]
87	REAL sulfate_pi(klon, klev) ! sulfate aerosol mass concentration [ug m-3] (pre-industrial value)
88	REAL sulfate(klon, klev) ! sulfate aerosol mass concentration [ug m-3
89	REAL cldtaupi(klon, klev) ! pre-industrial cloud opt thickness for diag
90	REAL re(klon, klev) ! cloud droplet effective radius [um]
91	REAL fl(klon, klev) ! xliq * rneb (denominator to re; fraction of liquid water clouds within the grid cell)
92	Real reliq(klon, klev), reice(klon, klev)
93
94
95	REAL tfi(klon,klev),qfi(klon,klev),t_seri(klon,klev)
96	REAL paprsfi(klon,klev+1),pplayfi(klon,klev)
97	integer read_climoz,read_oz
98	REAL wo(klon,klev),wo_dl(klon,klev)
99	REAL wo_ref(klon,klev),wo_dl_ref(klon,klev)
100	REAL wof1(klon,klev,1),wof2(klon,klev,2)
101	REAL wof1_ref(klon,klev,1),wof2_ref(klon,klev,2)
102	real dobson_u ! Dobson unit, in kg m-2
103	parameter (dobson_u=2.1415e-05)
104	real zmasse(klon, klev)
105	REAL albsfi(klon)
106	REAL tsolfi(klon)
107
108
109	Cjld integer l,ig0,m
110	integer l,ig0
111
112	REAL rlatfi(klon)
113
114	REAL cldliqfi(klon,klev),cldfrafi(klon,klev)
115	REAL cldtau(klon,klev), cldemi(klon,klev)
116
117	REAL heat(klon,klev),heat0(klon,klev),
118	$ cool(klon,klev),cool0(klon,klev)
119	REAL heat_ref(klon,klev),heat0_ref(klon,klev),
120	$ cool_ref(klon,klev),cool0_ref(klon,klev)
121	REAL d_heat_ini(klon,klev), d_heat0_ini(klon,klev),
122	$ d_cool_ini(klon,klev), d_cool0_ini(klon,klev)
123	REAL d_heat_adj(klon,klev), d_heat0_adj(klon,klev),
124	$ d_cool_adj(klon,klev), d_cool0_adj(klon,klev)
125
126	REAL topsw(klon),toplw(klon),solsw(klon),sollw(klon)
127	REAL topsw0(klon),toplw0(klon),solsw0(klon),sollw0(klon)
128	REAL sollwdown(klon),radsol(klon),albpla(klon)
129
130	REAL flx_lw_up(klon,klev+1),flx_lw_dn(klon,klev+1)
131	REAL flx_lw_up0(klon,klev+1),flx_lw_dn0(klon,klev+1) ! flux LW ciel clair up
132
133	REAL flx_sw_up(klon,klev+1),flx_sw_dn(klon,klev+1)
134	REAL flx_sw_up0(klon,klev+1),flx_sw_dn0(klon,klev+1) ! flux SW ciel clair up
135
136	REAL fdh(klon,klev) ! fixed dynamical heating
137	REAL dHrad_dT(klon,klev) ! derivative the radiative heating with temperature
138
139	c abderrahmane
140	c aerosols
141	real topswad(klon),solswad(klon),topswai(klon),solswai(klon) ! output: aerosol direct forcing at TOA and surface
142	real tau_ae(klon,klev,9,2),piz_ae(klon,klev,9,2)
143	real cg_ae(klon,klev,9,2) ! aerosol optical properties (see aeropt.F)
144	c ! (i.e., with a smaller droplet concentrationand thus larger droplet radii)
145	REAL bl95_b0, bl95_b1 ! Parameter in B&L 95-Formula
146	save bl95_b0, bl95_b1
147	logical ok_ade, ok_aie ! switches whether to use aerosol direct (indirect) effects or not
148	SAVE ok_ade, ok_aie
149
150	REAL aerindex(klon) ! POLDER aerosol index
151
152	LOGICAL new_aod
153	save new_aod
154	Real qsat(klon,klev),flwc(klon,klev), fiwc(klon,klev) ! Variable pour iflag_rrtm=1
155	Real topsw_aero(klon,9),topsw0_aero(klon,9)
156	Real solsw_aero(klon,9), solsw0_aero(klon,9)
157	Real topswad0_aero(klon), solswad0_aero(klon)
158	Real topswcf_aero(klon,3), solswcf_aero(klon,3)
159
160	logical ok_ade_ref, ok_aie_ref, ok_ade_modif, ok_aie_modif
161	save ok_ade_ref, ok_aie_ref, ok_ade_modif, ok_aie_modif
162
163
164	c reference profile
165	REAL tps_sw_ref(klon),tps_lw_ref(klon)
166	REAL tps_sw_ref0(klon),tps_lw_ref0(klon)
167	REAL toa_sw_ref(klon),toa_lw_ref(klon)
168	REAL toa_sw_ref0(klon),toa_lw_ref0(klon)
169	REAL srf_sw_ref(klon),srf_lw_ref(klon)
170	REAL srf_sw_ref0(klon),srf_lw_ref0(klon)
171	c
172	c initial values (with modified PARAMETER, but no strato. adjustment)
173	REAL tps_sw_ini(klon),tps_lw_ini(klon)
174	REAL tps_sw_ini0(klon),tps_lw_ini0(klon)
175	REAL toa_sw_ini(klon),toa_lw_ini(klon)
176	REAL toa_sw_ini0(klon),toa_lw_ini0(klon)
177	c initial anomalies
178	REAL d_tps_sw_ini(klon),d_tps_lw_ini(klon)
179	REAL d_tps_sw_ini0(klon),d_tps_lw_ini0(klon)
180	REAL d_toa_sw_ini(klon),d_toa_lw_ini(klon)
181	REAL d_toa_sw_ini0(klon),d_toa_lw_ini0(klon)
182	c
183	c adjusted values (with modified PARAMETER, and strato. adjustment)
184	REAL tps_sw_adj(klon),tps_lw_adj(klon)
185	REAL tps_sw_adj0(klon),tps_lw_adj0(klon)
186	REAL toa_sw_adj(klon),toa_lw_adj(klon)
187	REAL toa_sw_adj0(klon),toa_lw_adj0(klon)
188	REAL srf_sw_adj(klon),srf_lw_adj(klon)
189	REAL srf_sw_adj0(klon),srf_lw_adj0(klon)
190	c adjusted anomalies
191	REAL d_tps_sw_adj(klon),d_tps_lw_adj(klon)
192	REAL d_tps_sw_adj0(klon),d_tps_lw_adj0(klon)
193	REAL d_toa_sw_adj(klon),d_toa_lw_adj(klon)
194	REAL d_toa_sw_adj0(klon),d_toa_lw_adj0(klon)
195	REAL d_srf_sw_adj(klon),d_srf_lw_adj(klon)
196	REAL d_srf_sw_adj0(klon),d_srf_lw_adj0(klon)
197	c
198	c radiative budget of the stratosphere
199	REAL bilq_ref(klon),bilq_ini(klon),bilq_adj(klon)
200	!
201	REAL dT_tmp, dT_max(nbs),bilq_max(nbs)
202	c
203	INTEGER kpause(klon)
204	REAL ppause(klon)
205	c
206	REAL dist, rmu0(klon),fract(klon)
207	REAL zlongi
208	real solarlong0 ! si on veut imposer une valeur de la longitude solaire
209	data solarlong0/-999.99/
210	REAL x_ave, x_std, x_min, x_max
211	REAL undef
212	DATA undef/9999./
213
214	logical debut,ok_ozone
215	DATA debut/.true./
216	save debut
217
218	if (debut) THEN
219	c
220	print*,'Appel de suphel'
221	CALL suphel
222	!!
223	print*,'Lecture physiq.def'
224	!!
225	!Config Desc = Excentricite
226	!valeur AMIP II
227	R_ecc = 0.016715
228	call getin('R_ecc', R_ecc)
229	!Config Desc = Equinoxe
230	!valeur AMIP II
231	R_peri = 102.7
232	call getin('R_peri', R_peri)
233	!Config Desc = Inclinaison
234	!valeur AMIP II
235	R_incl = 23.441
236	call getin('R_incl', R_incl)
237	!Config Desc = Constante solaire en W/m2
238	!valeur AMIP II
239	solaire = 1365.
240	call getin('solaire', solaire)
241	!
242	! Proprietes des nuages
243	!
244	rad_froid = 35.0
245	call getin('rad_froid',rad_froid)
246	!
247	rad_chau1 = 13.0
248	call getin('rad_chau1',rad_chau1)
249	!
250	rad_chau2 = 9.0
251	call getin('rad_chau2',rad_chau2)
252	!
253	top_height = 3
254	call getin('top_height',top_height)
255	!
256	overlap = 3
257	call getin('overlap',overlap)
258	!
259	! Concentration des gaz a effet de serre
260	!
261	co2_ppm = 348.
262	call getin('co2_ppm', co2_ppm)
263	!
264	CH4_ppb = 1650.
265	call getin('CH4_ppb', CH4_ppb)
266	!
267	N2O_ppb=306.
268	call getin('N2O_ppb', N2O_ppb)
269	!
270	CFC11_ppt = 280.
271	call getin('CFC11_ppt',CFC11_ppt)
272	!
273	CFC12_ppt = 484.
274	call getin('CFC12_ppt',CFC12_ppt)
275	!
276	! Aerosols
277	!
278	!Config Desc = Aerosol direct effect not
279	ok_ade = .false.
280	! call getin('ok_ade', ok_ade)
281	!
282	!Config Desc = Aerosol indirect effect not
283	ok_aie = .false.
284	! call getin('ok_aie', ok_aie)
285	!
286	!Config Desc = Parameter in CDNC-maer link (Boucher&Lohmann 1995)
287	bl95_b0 = 1.7
288	call getin('bl95_b0', bl95_b0)
289	!
290	!Config Desc = Parameter in CDNC-maer link (Boucher&Lohmann 1995)
291	bl95_b1 = 0.2
292	call getin('bl95_b1', bl95_b1)
293	!
294	!Config Desc = Aerosol
295	! Test sur new_aod. Ce flag permet de retrouver les resultats de l'AR4
296	! il n'est utilisable que lors du couplage avec le SO4 seul
297	new_aod = .false.
298	call getin('new_aod', new_aod)
299	!
300	!read_climoz
301	read_climoz = 2
302	call getin('read_climoz', read_climoz)
303	if (read_climoz.ne.read_oz) then
304	print*,'Attention read_climoz de config.def
305	s different de celui correspondant a histmth.nc'
306	endif
307
308	t_glace_min = 258.
309	call getin('t_glace_min',t_glace_min)
310
311	t_glace_max = 273.13
312	call getin('t_glace_max',t_glace_max)
313
314	! les valeures lues sont prises comme valeures de références
315
316	solaire_ref = solaire
317	co2_ppm_ref = co2_ppm
318	CH4_ppb_ref = CH4_ppb
319	N2O_ppb_ref = N2O_ppb
320	CFC11_ppt_ref = CFC11_ppt
321	CFC12_ppt_ref = CFC12_ppt
322	ok_ade_ref = ok_ade
323	ok_aie_ref = ok_aie
324	!
325	RCO2_ref = co2_ppm_ref * 1.0e-06 * 44.011/28.97
326	RCH4_ref = CH4_ppb_ref * 1.0E-09 * 16.043/28.97
327	RN2O_ref = N2O_ppb_ref * 1.0E-09 * 44.013/28.97
328	RCFC11_ref = CFC11_ppt_ref * 1.0E-12 * 137.3686/28.97
329	RCFC12_ref = CFC12_ppt_ref * 1.0E-12 * 120.9140/28.97
330
331	PRINT*, "Valeurs de rad_chau1, rad_chau2, rad_froid :",
332	. rad_chau1, rad_chau2, rad_froid
333
334	PRINT*
335	PRINT*, "Lecture des valeures modifiées des paramètres"
336
337	read(,)solaire_modif
338	read(,)co2_ppm_modif
339	RCO2_modif = co2_ppm_modif * 1.0e-06 * 44.011/28.97
340
341	read(,)CH4_ppb_modif
342	RCH4_modif = CH4_ppb_modif * 1.0E-09 * 16.043/28.97
343
344	read(,)N2O_ppb_modif
345	RN2O_modif = N2O_ppb_modif * 1.0E-09 * 44.013/28.97
346
347	read(,)CFC11_ppt_modif
348	RCFC11_modif=CFC11_ppt_modif* 1.0E-12 * 137.3686/28.97
349
350	read(,)CFC12_ppt_modif
351	RCFC12_modif = CFC12_ppt_modif * 1.0E-12 * 120.9140/28.97
352
353	read(,)ok_ozone
354
355	ok_ade_modif=.false.
356	ok_aie_modif=.false.
357	!
358	write(,)' ##############################################'
359	write(,)' Valeures de références , perturbées et différences'
360	write(*,9000)' Constante solaire =',solaire_ref,solaire_modif
361	$ ,solaire_modif-solaire_ref
362	write(*,9000)' co2_ppm =',co2_ppm_ref,co2_ppm_modif
363	$ ,co2_ppm_modif-co2_ppm_ref
364	write(*,9000)' CH4_ppb =',CH4_ppb_ref,CH4_ppb_modif
365	$ ,CH4_ppb_modif-CH4_ppb_ref
366	write(*,9000)' N2O_ppb =',N2O_ppb_ref,N2O_ppb_modif
367	$ ,N2O_ppb_modif-N2O_ppb_ref
368	write(*,9000)' CFC11_ppt=',CFC11_ppt_ref,CFC11_ppt_modif
369	$ ,CFC11_ppt_modif-CFC11_ppt_ref
370	write(*,9000)' CFC12_ppt=',CFC12_ppt_ref,CFC12_ppt_modif
371	$ ,CFC12_ppt_modif-CFC12_ppt_ref
372	write(*,9001)' ok_ade=',ok_ade_ref, ok_ade_modif
373	write(*,9001)' ok_aie=',ok_aie_ref, ok_aie_modif
374	write(,)
375	9000 FORMAT (1x,A20,3(1pES15.6))
376	9001 FORMAT (1x,A20,2(L15))
377	! Initialisation rayonnrmrnt RRTM
378	! call iniradia(klon,klev,paprsfi(1,1:klev+1))
379	!
380	debut=.false.
381	endif
382	! calcul de zmasse
383	if (debug) print*,'rg=',rg
384	forall (k=1: klev) zmasse(:, k)=(paprsfi(:, k)-paprsfi(:, k+1))/rg
385	! Abderrahmane 26 11 2010
386	! A revfaire correctement
387	! Champ pour RRTM (a lire dans histmth ou recalculer)
388	! Pour l'instant on les met a 0
389	qsat=0.
390	flwc=0.
391	fiwc=0.
392	! initialiser 0 les sorties
393	! tropopause
394	tps_sw_ref=0.
395	tps_sw_ref0=0.
396	tps_lw_ref=0.
397	tps_lw_ref0=0.
398	tps_sw_ini=0.
399	tps_sw_ini0=0.
400	tps_lw_ini=0.
401	tps_lw_ini0=0.
402	d_tps_sw_ini=0.
403	d_tps_sw_ini0=0.
404	d_tps_lw_ini=0.
405	d_tps_lw_ini0=0.
406	d_tps_sw_adj=0.
407	d_tps_sw_adj0=0.
408	d_tps_lw_adj=0.
409	d_tps_lw_adj0=0.
410	! toa
411	toa_sw_ref=0.
412	toa_sw_ref0=0.
413	toa_lw_ref=0.
414	toa_lw_ref0=0.
415	toa_sw_ini=0.
416	toa_sw_ini0=0.
417	toa_lw_ini=0.
418	toa_lw_ini0=0.
419	d_toa_sw_ini=0.
420	d_toa_sw_ini0=0.
421	d_toa_lw_ini=0.
422	d_toa_lw_ini0=0.
423	d_toa_sw_adj=0.
424	d_toa_sw_adj0=0.
425	d_toa_lw_adj=0.
426	d_toa_lw_adj0=0.
427	! surface
428	srf_sw_ref=0.
429	srf_sw_ref0=0.
430	srf_lw_ref=0.
431	srf_lw_ref0=0.
432	d_srf_sw_adj=0.
433	d_srf_sw_adj0=0.
434	d_srf_lw_adj=0.
435	d_srf_lw_adj0=0.
436	!
437	! IF (debug) print*,'Appel ozonecm'
438	! wo(:, :, 1) = ozonecm(rlatfi, paprsfi, rjour)
439
440	IF (debug) print*,'Appel de tpause'
441	CALL tpause(tfi, pplayfi, kpause, ppause)
442
443	cCalculer epaisseur optique et emmissivite des nuages
444	IF (debug) print*,'Appel de newmicro'
445	CALL newmicro (paprsfi, pplayfi, ok_newmicro,
446	$ tfi, cldliqfi, cldfrafi, cldtau, cldemi,
447	. pch, pcl, pcm, pct, pctlwp,
448	s xflwp, xfiwp, xflwc, xfiwc,
449	e ok_aie,
450	e sulfate_ref, sulfate_pi,
451	e bl95_b0, bl95_b1,
452	s cldtaupi, re, fl, reliq, reice)
453
454
455	cpour un jour donne, calculer la longitude vraie de la terre
456	c (par rapport au point vernal-21 mars) dans son orbite solaire
457	c calculer aussi la distance terre-soleil (unite astronomique)
458	IF (debug) print*,'Appel de orbite, rjour=',rjour
459	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
460	! solarlong0
461	! !!! Atention, modif ad-hoc pour aquaplanette, quand on tourne un seul mois !!!!
462	! IF ( (m .eq. 1) .and. (rjour. eq.75.) ) solarlong0 = 0.
463	! IF (debug) print*,'m,rjour,solarlong0:',m,rjour,solarlong0
464	! if (solarlong0<-999.) then
465	CALL orbite(rjour,zlongi,dist)
466	if (debug) print*,'dist=',dist
467	! else
468	! zlongi=solarlong0 ! longitude solaire vraie
469	! dist=1. ! distance au soleil / moyenne
470	! write (,) 'Attention, orbite impose, solarlong0=',solarlong0
471	! endif
472	cCalculer la duree d'ensoleillement pour un jour et la hauteur
473	c du soleil (cosinus de l'angle zinithal) moyenne sur la journee
474	IF (debug) print*,'Appel angle'
475	CALL angle(zlongi,rlatfi,fract,rmu0)
476	IF (debug) print*,'zlongi:',zlongi,' rlatfi:',rlatfi
477	$ ,' fract:',fract,' rmu0:',rmu0
478	c
479
480	DO k = 1, klev
481	DO i=1,klon
482	t_seri(i,k) = tfi(i,k)
483	ENDDO
484	ENDDO
485	c
486	c Appel initial du rayonnement (reference)
487	! =======================================
488	c
489	IF (debug) print*,'Appel radlwsw no 1'
490
491	solaire=solaire_ref
492	RCO2=rco2_ref
493	RCH4=rch4_ref
494	RN2O=RN2O_ref
495	rcfc11=rcfc11_ref
496	rcfc12=rcfc12_ref
497	ok_ade=ok_ade_ref
498	ok_aie=ok_aie_ref
499
500	IF (debug) print*,'solaire, RCO2, ok_ade, ok_aie',
501	e solaire, RCO2, ok_ade, ok_aie
502	IF (debug) print*,'RCH4, RN2O, rcfc11, rcfc12 ',
503	e RCH4, RN2O, rcfc11, rcfc12
504	!
505	! Attention effet des aerosols dir et indir
506	! Abderrahmane ! A voir !
507	tau_ae=0.0
508	piz_ae=0.0
509	cg_ae=0.0
510
511	!Traitement cas ozone
512	if (read_climoz.eq.2) then
513	wof2_ref(:,:,1)=wo_ref(:,:)/dobson_u/1e3zmasse(:,:)rmo3/rmd
514	wof2_ref(:,:,2)=wo_dl_ref(:,:)/(dobson_u1e3/zmasse(:,:)/rmo3rmd)
515
516	if (debug) then
517	print*,'paprsfi(1,:)=',paprsfi(1,:)
518	print*,'pplayfi(1,:)=',pplayfi(1,:)
519	print*,'tsolfi(1)=',tsolfi(1)
520	print*,'albsfi(1)=',albsfi(1)
521	print*,'t_seri(1,:)=',t_seri(1,:)
522	print*,'qfi(1,:)=',qfi(1,:)
523	print*,'wof2_ref(1,:,:)=',wof2_ref(1,:,:)
524	print*,'cldfrafi(1,:)=',cldfrafi(1,:)
525	print*,'cldemi(1,:)=',cldemi(1,:)
526	print*,'cldtau(1,:)=',cldtau(1,:)
527	print*,'ok_ade, ok_aie=',ok_ade, ok_aie
528	endif
529
530	CALL radlwsw(dist, rmu0, fract,
531	e paprsfi, pplayfi,tsolfi,albsfi,
532	e albsfi,t_seri,qfi,wof2_ref,
533	e cldfrafi, cldemi, cldtau,
534	e ok_ade, ok_aie,
535	e tau_ae, piz_ae, cg_ae,
536	e cldtaupi, new_aod,
537	N qsat, flwc, fiwc,
538	s heat,heat0,cool,cool0,radsol,albpla,
539	s topsw,toplw,solsw,sollw,
540	a sollwdown,
541	s topsw0,toplw0,solsw0,sollw0,
542	s flx_lw_dn0, flx_lw_dn, flx_lw_up0, flx_lw_up,
543	a flx_sw_dn0, flx_sw_dn, flx_sw_up0, flx_sw_up,
544	a topswad, solswad,
545	a topswai, solswai,
546	a topswad0_aero, solswad0_aero,
547	N topsw_aero, topsw0_aero,
548	N solsw_aero, solsw0_aero,
549	N topswcf_aero, solswcf_aero )
550	else
551	wof1_ref(:,:,1)=wo_ref(:,:)/(dobson_u1e3/zmasse(:,:)/rmo3rmd)
552	CALL radlwsw(dist, rmu0, fract,
553	e paprsfi, pplayfi,tsolfi,albsfi,
554	e albsfi,t_seri,qfi,wof1_ref,
555	e cldfrafi, cldemi, cldtau,
556	e ok_ade, ok_aie,
557	e tau_ae, piz_ae, cg_ae,
558	e cldtaupi, new_aod,
559	N qsat, flwc, fiwc,
560	s heat,heat0,cool,cool0,radsol,albpla,
561	s topsw,toplw,solsw,sollw,
562	a sollwdown,
563	s topsw0,toplw0,solsw0,sollw0,
564	s flx_lw_dn0, flx_lw_dn, flx_lw_up0, flx_lw_up,
565	a flx_sw_dn0, flx_sw_dn, flx_sw_up0, flx_sw_up,
566	a topswad, solswad,
567	a topswai, solswai,
568	a topswad0_aero, solswad0_aero,
569	N topsw_aero, topsw0_aero,
570	N solsw_aero, solsw0_aero,
571	N topswcf_aero, solswcf_aero )
572	endif
573
574	c
575	DO k = 1, klev !diagnostiquer le FDH
576	DO i = 1,klon
577	fdh(i,k) = -(heat(i,k)-cool(i,k))
578	ENDDO
579	ENDDO
580	!
581	! Flux net
582	DO i = 1, klon
583	! Ciel avec nuages
584	! SW
585	tps_sw_ref(i)=flx_sw_dn(i,kpause(i))-flx_sw_up(i,kpause(i))
586	toa_sw_ref(i)=flx_sw_dn(i,klev+1)-flx_sw_up(i,klev+1)
587	srf_sw_ref(i)=flx_sw_dn(i,1)-flx_sw_up(i,1)
588	! LW
589	tps_lw_ref(i)=-1.*(flx_lw_dn(i,kpause(i))+flx_lw_up(i,kpause(i)))
590	toa_lw_ref(i)=-1.*(flx_lw_dn(i,klev+1)+flx_lw_up(i,klev+1))
591	srf_lw_ref(i)=-1.*(flx_lw_dn(i,1)+flx_lw_up(i,1))
592	!
593	! ciel clair
594	! SW
595	tps_sw_ref0(i)=flx_sw_dn0(i,kpause(i))-flx_sw_up0(i,kpause(i))
596	toa_sw_ref0(i)=flx_sw_dn0(i,klev+1)-flx_sw_up0(i,klev+1)
597	srf_sw_ref0(i)=flx_sw_dn0(i,1)-flx_sw_up0(i,1)
598	! LW
599	tps_lw_ref0(i)=-1.*(flx_lw_dn0(i,kpause(i))
600	$ +flx_lw_up0(i,kpause(i)))
601	toa_lw_ref0(i)=-1.*(flx_lw_dn0(i,klev+1)+flx_lw_up0(i,klev+1))
602	srf_lw_ref0(i)=-1.*(flx_lw_dn0(i,1)+flx_lw_up0(i,1))
603	ENDDO
604	!
605	! reference radiative inbalance of the stratosphere
606	DO i = 1, klon
607	bilq_ref(i) = (toa_sw_ref(i)-tps_sw_ref(i))
608	$ +(toa_lw_ref(i)-tps_lw_ref(i))
609	END DO
610	! Bilan rad dans l'atmosphere:
611	heat_ref(:,:) = heat(:,:)
612	cool_ref(:,:) = cool(:,:)
613	heat0_ref(:,:) = heat0(:,:)
614	cool0_ref(:,:) = cool0(:,:)
615	!
616	! End of computation of the reference case
617	!
618	! Calcul de la derivee du bilan en fonction de la temperature
619	! ===========================================================
620	!
621	! On incremente la temp de 1 dans la tropopause
622	DO k = 1, klev
623	DO i = 1, klon
624	IF (k.GE.kpause(i)) t_seri(i,k)=t_seri(i,k)+1.
625	ENDDO
626	ENDDO
627
628	if (read_climoz.eq.2) then
629	wof2_ref(:,:,1)=wo_ref(:,:)/(dobson_u1e3/zmasse(:,:)/rmo3rmd)
630	wof2_ref(:,:,2)=wo_dl_ref(:,:)/(dobson_u1e3/zmasse(:,:)/rmo3rmd)
631	CALL radlwsw(dist, rmu0, fract,
632	e paprsfi, pplayfi,tsolfi,albsfi,
633	e albsfi,t_seri,qfi,wof2_ref,
634	e cldfrafi, cldemi, cldtau,
635	e ok_ade, ok_aie,
636	e tau_ae, piz_ae, cg_ae,
637	e cldtaupi, new_aod,
638	N qsat, flwc, fiwc,
639	s heat,heat0,cool,cool0,radsol,albpla,
640	s topsw,toplw,solsw,sollw,
641	a sollwdown,
642	s topsw0,toplw0,solsw0,sollw0,
643	s flx_lw_dn0, flx_lw_dn, flx_lw_up0, flx_lw_up,
644	a flx_sw_dn0, flx_sw_dn, flx_sw_up0, flx_sw_up,
645	a topswad, solswad,
646	a topswai, solswai,
647	a topswad0_aero, solswad0_aero,
648	N topsw_aero, topsw0_aero,
649	N solsw_aero, solsw0_aero,
650	N topswcf_aero, solswcf_aero )
651	else
652	wof1_ref(:,:,1)=wo_ref(:,:)/(dobson_u1e3/zmasse(:,:)/rmo3rmd)
653	CALL radlwsw(dist, rmu0, fract,
654	e paprsfi, pplayfi,tsolfi,albsfi,
655	e albsfi,t_seri,qfi,wof1_ref,
656	e cldfrafi, cldemi, cldtau,
657	e ok_ade, ok_aie,
658	e tau_ae, piz_ae, cg_ae,
659	e cldtaupi, new_aod,
660	N qsat, flwc, fiwc,
661	s heat,heat0,cool,cool0,radsol,albpla,
662	s topsw,toplw,solsw,sollw,
663	a sollwdown,
664	s topsw0,toplw0,solsw0,sollw0,
665	s flx_lw_dn0, flx_lw_dn, flx_lw_up0, flx_lw_up,
666	a flx_sw_dn0, flx_sw_dn, flx_sw_up0, flx_sw_up,
667	a topswad, solswad,
668	a topswai, solswai,
669	a topswad0_aero, solswad0_aero,
670	N topsw_aero, topsw0_aero,
671	N solsw_aero, solsw0_aero,
672	N topswcf_aero, solswcf_aero )
673	endif
674	!
675	dHrad_dT(:,:)=-1.
676	DO k = 1, klev !diagnostiquer le d FDH / d T
677	DO i = 1,klon
678	IF (k.GE.kpause(i))
679	$ dHrad_dT(i,k) = (heat(i,k)-cool(i,k)) + fdh(i,k)
680	ENDDO
681	ENDDO
682	! On remet la valeure initial du profil de temperature
683	DO k = 1, klev
684	DO i=1,klon
685	t_seri(i,k) = tfi(i,k)
686	ENDDO
687	ENDDO
688	! Fin du calcul de la derivee du bilan en fonction de la temperature
689	!
690	! Calcul des flux radiatifs apres perturbation
691	! ============================================
692	!
693	solaire=solaire_modif
694	RCO2=rco2_modif
695	RCH4=rch4_modif
696	RN2O=RN2O_modif
697	rcfc11=rcfc11_modif
698	rcfc12=rcfc12_modif
699	ok_ade=ok_ade_modif
700	ok_aie=ok_aie_modif
701
702	if(.NOT.ok_ozone)then
703	wo=wo_ref
704	wo_dl=wo_dl_ref
705	endif
706
707	!
708	! On itere pour ajuster la stratosphere
709	DO 88888 itap = 1, nbs
710
711	IF (debug) print*,'Appel radlwsw no 2'
712	IF (debug) print*,'solaire, RCO2, ok_ade, ok_aie',
713	e solaire, RCO2, ok_ade, ok_aie
714
715	if (read_climoz.eq.2) then
716	wof2(:,:,1)=wo(:,:)/(dobson_u1e3/zmasse(:,:)/rmo3rmd)
717	wof2(:,:,2)=wo_dl(:,:)/(dobson_u1e3/zmasse(:,:)/rmo3rmd)
718	CALL radlwsw(dist, rmu0, fract,
719	e paprsfi, pplayfi,tsolfi,albsfi,
720	e albsfi,t_seri,qfi,wof2,
721	e cldfrafi, cldemi, cldtau,
722	e ok_ade, ok_aie,
723	e tau_ae, piz_ae, cg_ae,
724	e cldtaupi, new_aod,
725	N qsat, flwc, fiwc,
726	s heat,heat0,cool,cool0,radsol,albpla,
727	s topsw,toplw,solsw,sollw,
728	a sollwdown,
729	s topsw0,toplw0,solsw0,sollw0,
730	s flx_lw_dn0, flx_lw_dn, flx_lw_up0, flx_lw_up,
731	a flx_sw_dn0, flx_sw_dn, flx_sw_up0, flx_sw_up,
732	a topswad, solswad,
733	a topswai, solswai,
734	a topswad0_aero, solswad0_aero,
735	N topsw_aero, topsw0_aero,
736	N solsw_aero, solsw0_aero,
737	N topswcf_aero, solswcf_aero )
738	else
739	wof1(:,:,1)=wo(:,:)/(dobson_u1e3/zmasse(:,:)/rmo3rmd)
740	CALL radlwsw(dist, rmu0, fract,
741	e paprsfi, pplayfi,tsolfi,albsfi,
742	e albsfi,t_seri,qfi,wof1,
743	e cldfrafi, cldemi, cldtau,
744	e ok_ade, ok_aie,
745	e tau_ae, piz_ae, cg_ae,
746	e cldtaupi, new_aod,
747	N qsat, flwc, fiwc,
748	s heat,heat0,cool,cool0,radsol,albpla,
749	s topsw,toplw,solsw,sollw,
750	a sollwdown,
751	s topsw0,toplw0,solsw0,sollw0,
752	s flx_lw_dn0, flx_lw_dn, flx_lw_up0, flx_lw_up,
753	a flx_sw_dn0, flx_sw_dn, flx_sw_up0, flx_sw_up,
754	a topswad, solswad,
755	a topswai, solswai,
756	a topswad0_aero, solswad0_aero,
757	N topsw_aero, topsw0_aero,
758	N solsw_aero, solsw0_aero,
759	N topswcf_aero, solswcf_aero )
760	endif
761	c Ajouter la tendance des rayonnements et le FDH
762	c
763	dT_max(itap) = 0.
764	DO k = 1, klev
765	DO i = 1, klon
766	IF (k.GE.kpause(i)) THEN
767	dT_tmp = (heat(i,k)-cool(i,k)+fdh(i,k))
768	$ /(86400./dtime-0.5*dHrad_dT(i,k))
769	! print *,i,k, dT_tmp,dHrad_dT(i,k)
770	! dT_tmp = (heat(i,k)-cool(i,k)+fdh(i,k))
771	t_seri(i,k) = t_seri(i,k) + dT_tmp
772	dT_max(itap) = max(dT_max(itap),abs(dT_tmp))
773	ENDIF
774	ENDDO
775	ENDDO
776	c
777	IF ( itap .EQ. 1 ) THEN
778	DO i = 1, klon
779	tps_sw_ini(i) = flx_sw_dn(i,kpause(i))
780	. - flx_sw_up(i,kpause(i))
781	toa_sw_ini(i) = flx_sw_dn(i,klev+1)
782	. - flx_sw_up(i,klev+1)
783
784	tps_lw_ini(i) = -1.*(flx_lw_dn(i,kpause(i))
785	. + flx_lw_up(i,kpause(i)))
786	toa_lw_ini(i) = -1.*(flx_lw_dn(i,klev+1)
787	. + flx_lw_up(i,klev+1))
788
789	c difference relative to the reference values
790	d_tps_sw_ini(i) = tps_sw_ini(i)-tps_sw_ref(i)
791	d_toa_sw_ini(i) = toa_sw_ini(i)-toa_sw_ref(i)
792	d_tps_lw_ini(i) = tps_lw_ini(i)-tps_lw_ref(i)
793	d_toa_lw_ini(i) = toa_lw_ini(i)-toa_lw_ref(i)
794	C
795	tps_sw_ini0(i) = flx_sw_dn0(i,kpause(i))
796	. - flx_sw_up0(i,kpause(i))
797	toa_sw_ini0(i) = flx_sw_dn0(i,klev+1)
798	. - flx_sw_up0(i,klev+1)
799
800	tps_lw_ini0(i) = -1.*(flx_lw_dn0(i,kpause(i))
801	. + flx_lw_up0(i,kpause(i)))
802	toa_lw_ini0(i) = -1.*(flx_lw_dn0(i,klev+1)
803	. + flx_lw_up0(i,klev+1))
804	c difference relative to the reference values
805	d_tps_sw_ini0(i) = tps_sw_ini0(i)-tps_sw_ref0(i)
806	d_toa_sw_ini0(i) = toa_sw_ini0(i)-toa_sw_ref0(i)
807	d_tps_lw_ini0(i) = tps_lw_ini0(i)-tps_lw_ref0(i)
808	d_toa_lw_ini0(i) = toa_lw_ini0(i)-toa_lw_ref0(i)
809	END DO
810	!
811	IF (debug) then
812	CALL cstat(klon,toa_sw_ini,x_ave,x_std,x_min,x_max,undef)
813	WRITE (*,9002) "toa_sw_ini",x_ave,x_std,x_min,x_max
814	CALL cstat(klon,tps_sw_ini,x_ave,x_std,x_min,x_max,undef)
815	WRITE (*,9002) "tps_sw_ini",x_ave,x_std,x_min,x_max
816	CALL cstat(klon,toa_sw_ini0,x_ave,x_std,x_min,x_max,undef)
817	WRITE (*,9002) "toa_sw_ini0",x_ave,x_std,x_min,x_max
818	CALL cstat(klon,tps_sw_ini0,x_ave,x_std,x_min,x_max,undef)
819	WRITE (*,9002) "tps_sw_ini0",x_ave,x_std,x_min,x_max
820	!
821	CALL cstat(klon,toa_lw_ini,x_ave,x_std,x_min,x_max,undef)
822	WRITE (*,9002) "toa_lw_ini",x_ave,x_std,x_min,x_max
823	CALL cstat(klon,tps_lw_ini,x_ave,x_std,x_min,x_max,undef)
824	WRITE (*,9002) "tps_lw_ini",x_ave,x_std,x_min,x_max
825	CALL cstat(klon,toa_lw_ini0,x_ave,x_std,x_min,x_max,undef)
826	WRITE (*,9002) "toa_lw_ini0",x_ave,x_std,x_min,x_max
827	CALL cstat(klon,tps_lw_ini0,x_ave,x_std,x_min,x_max,undef)
828	WRITE (*,9002) "tps_lw_ini0",x_ave,x_std,x_min,x_max
829	!
830	END IF
831	9002 FORMAT (1x,A12,10(1pE13.6))
832	!
833	! initial radiatiave inbalance of the stratosphere
834	DO i = 1, klon
835	bilq_ini(i) = (toa_sw_ini(i)-tps_sw_ini(i))
836	$ + (toa_lw_ini(i)-tps_lw_ini(i))
837	END DO
838	!
839	! Variation initiale du bilan rad dans l'atmosphere:
840	d_heat_ini(:,:) = heat(:,:) - heat_ref(:,:)
841	d_cool_ini(:,:) = cool(:,:) - cool_ref(:,:)
842	d_heat0_ini(:,:) = heat0(:,:) - heat0_ref(:,:)
843	d_cool0_ini(:,:) = cool0(:,:) - cool0_ref(:,:)
844	!
845	END IF ! ( itap .EQ. 1 )
846	c
847	88888 CONTINUE
848
849	DO i = 1, klon
850	!
851	tps_sw_adj(i) = flx_sw_dn(i,kpause(i))
852	. - flx_sw_up(i,kpause(i))
853	toa_sw_adj(i) = flx_sw_dn(i,klev+1)
854	. - flx_sw_up(i,klev+1)
855	srf_sw_adj(i) = flx_sw_dn(i,1)
856	. - flx_sw_up(i,1)
857	tps_lw_adj(i) = -1.*(flx_lw_dn(i,kpause(i))
858	. + flx_lw_up(i,kpause(i)))
859	toa_lw_adj(i) = -1.*(flx_lw_dn(i,klev+1)
860	. + flx_lw_up(i,klev+1))
861	srf_lw_adj(i) = -1.*(flx_lw_dn(i,1)
862	. + flx_lw_up(i,1))
863	c difference relative to the reference values
864	d_tps_sw_adj(i) = tps_sw_adj(i)-tps_sw_ref(i)
865	d_toa_sw_adj(i) = toa_sw_adj(i)-toa_sw_ref(i)
866	d_srf_sw_adj(i) = srf_sw_adj(i)-srf_sw_ref(i)
867	d_tps_lw_adj(i) = tps_lw_adj(i)-tps_lw_ref(i)
868	d_toa_lw_adj(i) = toa_lw_adj(i)-toa_lw_ref(i)
869	d_srf_lw_adj(i) = srf_lw_adj(i)-srf_lw_ref(i)
870	C
871	tps_sw_adj0(i) = flx_sw_dn0(i,kpause(i))
872	. - flx_sw_up0(i,kpause(i))
873	toa_sw_adj0(i) = flx_sw_dn0(i,klev+1)
874	. - flx_sw_up0(i,klev+1)
875	srf_sw_adj0(i) = flx_sw_dn0(i,1)
876	. - flx_sw_up0(i,1)
877	tps_lw_adj0(i) = -1.*(flx_lw_dn0(i,kpause(i))
878	. + flx_lw_up0(i,kpause(i)))
879	toa_lw_adj0(i) = -1.*(flx_lw_dn0(i,klev+1)
880	. + flx_lw_up0(i,klev+1))
881	srf_lw_adj0(i) = -1.*(flx_lw_dn0(i,1)
882	. + flx_lw_up0(i,1))
883	c difference relative to the reference values
884	d_tps_sw_adj0(i) = tps_sw_adj0(i)-tps_sw_ref0(i)
885	d_toa_sw_adj0(i) = toa_sw_adj0(i)-toa_sw_ref0(i)
886	d_srf_sw_adj0(i) = srf_sw_adj0(i)-srf_sw_ref0(i)
887	d_tps_lw_adj0(i) = tps_lw_adj0(i)-tps_lw_ref0(i)
888	d_toa_lw_adj0(i) = toa_lw_adj0(i)-toa_lw_ref0(i)
889	d_srf_lw_adj0(i) = srf_lw_adj0(i)-srf_lw_ref0(i)
890	ENDDO
891	c adjusted radiatiave inbalance of the stratosphere
892	DO i = 1, klon
893	bilq_adj(i) = (toa_sw_adj(i)-tps_sw_adj(i))
894	$ + (toa_lw_adj(i)-tps_lw_adj(i))
895	END DO
896	!
897	! Variation ajustee du bilan rad dans l'atmosphere:
898	d_heat_adj(:,:) = heat(:,:) - heat_ref(:,:)
899	d_cool_adj(:,:) = cool(:,:) - cool_ref(:,:)
900	d_heat0_adj(:,:) = heat0(:,:) - heat0_ref(:,:)
901	d_cool0_adj(:,:) = cool0(:,:) - cool0_ref(:,:)
902	C
903	WRITE (*,9010) 'dT_max',(dT_max(itap),itap=1,nbs)
904	9010 FORMAT (1x,A8,(10E13.6))
905	c
906	RETURN
907	END
908	c======================================================================
909	C
910	SUBROUTINE tpause(t, pls, kpause, ppause)
911	use dimphy
912	IMPLICIT none
913	c
914	REAL plsmin
915	PARAMETER (plsmin=35000.0)
916	REAL seuil
917	PARAMETER (seuil=-2.0e-3)
918	ccc PARAMETER (seuil=2.0e-3)
919	c
920	#include "dimensions.h"
921
922	REAL t(klon,klev)
923	REAL pls(klon,klev)
924	INTEGER kpause(klon)
925	REAL ppause(klon)
926	c
927	LOGICAL deja(klon)
928	INTEGER kmin(klon)
929	INTEGER i, k
930	REAL deltat, deltaz
931
932	c
933	DO i = 1, klon
934	deja(i) = .FALSE.
935	ENDDO
936	c
937	DO k = 1, klev
938	DO i = 1, klon
939	IF (.NOT.deja(i) .AND. pls(i,k).LE.plsmin) THEN
940	kmin(i) = k
941	deja(i) = .TRUE.
942	ENDIF
943	ENDDO
944	ENDDO
945	c
946	DO i = 1, klon
947	deja(i) = .FALSE.
948	kpause(i) = kmin(i)
949	ppause(i) = (pls(i,kmin(i))+pls(i,kmin(i)-1))/2.0
950	ENDDO
951	c
952	DO k = 1, klev
953	DO i = 1, klon
954	IF (.NOT.deja(i) .AND. k.GE.kmin(i)) THEN
955	deltat = t(i,k) - t(i,k-1)
956	deltaz = 287.0/9.81*(t(i,k)+t(i,k-1))/(pls(i,k)+pls(i,k-1))
957	. *(pls(i,k-1)-pls(i,k))
958	IF (deltat/deltaz .GT. seuil) THEN
959	C JLD, 2009/11/03: on dit que la limite se situe a la maille du dessous
960	C kpause(i) = k
961	kpause(i) = k-1
962	CJLD ppause(i) = (pls(i,k)+pls(i,k-1))/2.0
963	ppause(i) = pls(i,k-1)
964	deja(i) = .TRUE.
965	ENDIF
966	ENDIF
967	ENDDO
968	ENDDO
969	c
970	RETURN
971	END
972	C
973	c======================================================================
974	C
975	SUBROUTINE lirehist(tapename,iim,jjmp1,klev,ngridmx,mois,m,reado,
976	. tsol,paprs,pplay,t,q,cldfra,cldwcon,ozone,ozone_dl,
977	. rlon, rlat, phis, aire, albs,
978	. rlon_1d, rlat_1d, presnivs)
979	IMPLICIT none
980	c
981	CHARACTER() tapename
982	INTEGER iim, jjmp1, klev, mois, m
983	INTEGER ngridmx
984	REAL rlat(ngridmx), rlon(ngridmx), presnivs(klev)
985	REAL phis(ngridmx), aire(ngridmx), albs(ngridmx)
986	REAL tsol(ngridmx)
987	REAL paprs(ngridmx,klev+1)
988	REAL pplay(ngridmx,klev), t(ngridmx,klev), q(ngridmx,klev)
989	REAL cldfra(ngridmx,klev)
990	REAL cldwcon(ngridmx,klev) ! cloud water content (ice+liquide)
991	integer reado
992	REAL ozone(ngridmx,klev),ozone_dl(ngridmx,klev)
993	REAL cldliq(ngridmx,klev),cldice(ngridmx,klev)
994	REAL psol(ngridmx)
995	C
996	LOGICAL debug
997	PARAMETER (debug=.false.)
998	C PARAMETER (debug=.true.)
999	REAL solsdn(ngridmx), solsup(ngridmx)
1000	REAL rlat_1d(jjmp1), rlon_1d(iim)
1001	c
1002	INTEGER debut4(4), epais4(4), debut3(3), epais3(3)
1003	C
1004	REAL zx_tmp_2d(iim,jjmp1)
1005	REAL zx_tmp_3d(iim,jjmp1,klev)
1006	c
1007	#include "netcdf.inc"
1008	c
1009	INTEGER nid, varid, ierr
1010	INTEGER i, j, k, l, n
1011	c
1012	IF (debug) PRINT *,'Entree de lirehist'
1013	ierr = NF_OPEN (tapename, NF_NOWRITE,nid)
1014	IF (ierr.NE.NF_NOERR) THEN
1015	PRINT*, "Probleme d ouverture du fichier:"//tapename
1016	CALL abort
1017	ENDIF
1018	c
1019	c La taille du fichier en longitude:
1020	c
1021	IF (debug) PRINT *,'lecture lon'
1022	ierr = NF_INQ_DIMID(nid,"lon",varid)
1023	IF (ierr.NE.NF_NOERR) THEN
1024	PRINT*, "Probleme pour avoir identificateur de lon"
1025	CALL abort
1026	ENDIF
1027	ierr = NF_INQ_DIMLEN(nid,varid,i)
1028	IF (ierr.NE.NF_NOERR) THEN
1029	PRINT*, "Probleme pour avoir la taille de lon"
1030	CALL abort
1031	ENDIF
1032	IF (i.NE.iim) THEN
1033	PRINT*, "Dimension de lon fausse:", i, iim
1034	CALL abort
1035	ENDIF
1036	c
1037	c La taille du fichier en latitude:
1038	c
1039	IF (debug) PRINT *,'lecture lat'
1040	ierr = NF_INQ_DIMID(nid,"lat",varid)
1041	IF (ierr.NE.NF_NOERR) THEN
1042	PRINT*, "Probleme pour avoir identificateur de lat"
1043	CALL abort
1044	ENDIF
1045	ierr = NF_INQ_DIMLEN(nid,varid,j)
1046	IF (ierr.NE.NF_NOERR) THEN
1047	PRINT*, "Probleme pour avoir la taille de lat"
1048	CALL abort
1049	ENDIF
1050	IF (j.NE.jjmp1) THEN
1051	PRINT*, "Dimension lat fausse:", j, jjmp1
1052	CALL abort
1053	ENDIF
1054	c
1055	c La taille du fichier en vertical:
1056	c
1057	IF (debug) PRINT *,'lecture dimension presniv'
1058	ierr = NF_INQ_DIMID(nid,"presnivs",varid)
1059	IF (ierr.NE.NF_NOERR) THEN
1060	PRINT*, "Probleme pour avoir identificateur de z"
1061	CALL abort
1062	ENDIF
1063	ierr = NF_INQ_DIMLEN(nid,varid,k)
1064	IF (ierr.NE.NF_NOERR) THEN
1065	PRINT*, "Probleme pour avoir la taille de z"
1066	CALL abort
1067	ENDIF
1068	IF (k.NE.klev) THEN
1069	PRINT*, "Dimension z fausse:", k, klev
1070	CALL abort
1071	ENDIF
1072	c
1073	c La taille du fichier en temps:
1074	c
1075	IF (debug) PRINT *,'lecture dimension time'
1076	ierr = NF_INQ_DIMID(nid,"time_counter",varid)
1077	IF (ierr.NE.NF_NOERR) THEN
1078	PRINT*, "Probleme pour avoir identificateur de time_counter"
1079	CALL abort
1080	ENDIF
1081	ierr = NF_INQ_DIMLEN(nid,varid,l)
1082	IF (ierr.NE.NF_NOERR) THEN
1083	PRINT*, "Probleme pour avoir la taille de time_counter"
1084	CALL abort
1085	ENDIF
1086	IF (m.GT.l) THEN
1087	PRINT*, "Dimension time_counter fausse:", l, m
1088	CALL abort
1089	ENDIF
1090	IF (m.GT.mois) THEN
1091	PRINT*, "Numero mois sup. dim des mois:", m, mois
1092	CALL abort
1093	ENDIF
1094	c
1095	IF (debug) PRINT *,'lecture presniv'
1096	ierr = NF_INQ_VARID (nid,"presnivs",varid)
1097	IF (ierr.NE.NF_NOERR) THEN
1098	PRINT*, "Probleme pour avoir identificateur presnivs"
1099	CALL abort
1100	ENDIF
1101
1102	#ifdef NC_DOUBLE
1103	ierr = NF_GET_VAR_DOUBLE(nid, varid, presnivs)
1104	#else
1105	ierr = NF_GET_VAR_REAL(nid, varid, presnivs)
1106	#endif
1107	IF (ierr.NE.NF_NOERR) THEN
1108	PRINT*, "Probleme pour extraire presnivs"
1109	CALL abort
1110	ENDIF
1111	c ----
1112	IF (debug) PRINT *,'lecture lat'
1113	ierr = NF_INQ_VARID (nid,"lat",varid)
1114	IF (ierr.NE.NF_NOERR) THEN
1115	PRINT*, "Probleme pour avoir identificateur lat"
1116	CALL abort
1117	ENDIF
1118
1119	#ifdef NC_DOUBLE
1120	ierr = NF_GET_VAR_DOUBLE(nid, varid, rlat_1d)
1121	#else
1122	ierr = NF_GET_VAR_REAL(nid, varid, rlat_1d)
1123	#endif
1124
1125	IF (ierr.NE.NF_NOERR) THEN
1126	PRINT*, "Probleme pour extraire lat"
1127	CALL abort
1128	ENDIF
1129	DO i=1,iim
1130	zx_tmp_2d(i,:)=rlat_1d(:)
1131	END DO
1132	CALL gr_fi_lu(1,ngridmx,iim,jjmp1,rlat,zx_tmp_2d)
1133	c ----
1134	IF (debug) PRINT *,'lecture lon'
1135	ierr = NF_INQ_VARID (nid,"lon",varid)
1136	IF (ierr.NE.NF_NOERR) THEN
1137	PRINT*, "Probleme pour avoir identificateur lon"
1138	CALL abort
1139	ENDIF
1140	#ifdef NC_DOUBLE
1141	ierr = NF_GET_VAR_DOUBLE(nid, varid, rlon_1d)
1142	#else
1143	ierr = NF_GET_VAR_REAL(nid, varid, rlon_1d)
1144	#endif
1145
1146	IF (ierr.NE.NF_NOERR) THEN
1147	PRINT*, "Probleme pour extraire lon"
1148	CALL abort
1149	ENDIF
1150	DO j=1,jjmp1
1151	zx_tmp_2d(:,j)=rlon_1d(:)
1152	END DO
1153	CALL gr_fi_lu(1,ngridmx,iim,jjmp1,rlon,zx_tmp_2d)
1154	c
1155	c Dans certains fichier hist, phis et aire ont une DIMENSION du temps
1156	debut3(1) = 1
1157	debut3(2) = 1
1158	debut3(3) =1
1159	epais3(1) = iim
1160	epais3(2) = jjmp1
1161	epais3(3) = 1
1162	C ----
1163	IF (debug) PRINT *,'lecture phis'
1164	ierr = NF_INQ_VARID (nid,"phis",varid)
1165	IF (ierr.NE.NF_NOERR) THEN
1166	PRINT*, "Probleme pour avoir identificateur phis"
1167	CALL abort
1168	ENDIF
1169	#ifdef NC_DOUBLE
1170	CJLD ierr = NF_GET_VAR_DOUBLE(nid, varid, phis)
1171	ierr = NF_GET_VARA_DOUBLE(nid,varid,debut3,epais3,zx_tmp_2d)
1172	#else
1173	CJLD ierr = NF_GET_VAR_REAL(nid, varid, phis)
1174	ierr = NF_GET_VARA_REAL(nid,varid,debut3,epais3,zx_tmp_2d)
1175	#endif
1176
1177	IF (ierr.NE.NF_NOERR) THEN
1178	PRINT*, "Probleme pour extraire phis"
1179	CALL abort
1180	ENDIF
1181	CALL gr_fi_lu(1,ngridmx,iim,jjmp1,phis,zx_tmp_2d)
1182	C
1183	IF (debug) PRINT *,'lecture aire'
1184	ierr = NF_INQ_VARID (nid,"aire",varid)
1185	IF (ierr.NE.NF_NOERR) THEN
1186	PRINT*, "Probleme pour avoir identificateur aire"
1187	CALL abort
1188	ENDIF
1189	#ifdef NC_DOUBLE
1190	CJLD ierr = NF_GET_VAR_DOUBLE(nid, varid, aire)
1191	ierr = NF_GET_VARA_DOUBLE(nid,varid,debut3,epais3,zx_tmp_2d)
1192	#else
1193	CJLD ierr = NF_GET_VAR_REAL(nid, varid, aire)
1194	ierr = NF_GET_VARA_REAL(nid,varid,debut3,epais3,zx_tmp_2d)
1195	#endif
1196
1197	IF (ierr.NE.NF_NOERR) THEN
1198	PRINT*, "Probleme pour extraire aire"
1199	CALL abort
1200	ENDIF
1201	CALL gr_fi_lu(1,ngridmx,iim,jjmp1,aire,zx_tmp_2d)
1202	c ----
1203	debut3(1) = 1
1204	debut3(2) = 1
1205	debut3(3) = m
1206	epais3(1) = iim
1207	epais3(2) = jjmp1
1208	epais3(3) = 1
1209	c
1210	IF (debug) PRINT *,'lecture tsol'
1211	ierr = NF_INQ_VARID(nid,"tsol", varid)
1212	IF (ierr.NE.NF_NOERR) THEN
1213	PRINT*, "Probleme pour avoir identificateur de tsol"
1214	CALL abort
1215	ENDIF
1216	#ifdef NC_DOUBLE
1217	ierr = NF_GET_VARA_DOUBLE(nid,varid,debut3,epais3,zx_tmp_2d)
1218	#else
1219	ierr = NF_GET_VARA_REAL(nid,varid,debut3,epais3,zx_tmp_2d)
1220	#endif
1221	IF (ierr.NE.NF_NOERR) THEN
1222	PRINT*, "Probleme pour extraire tsol"
1223	CALL abort
1224	ENDIF
1225	CALL gr_fi_lu(1,ngridmx,iim,jjmp1,tsol,zx_tmp_2d)
1226	c ----
1227	IF (debug) PRINT *,'lecture psol'
1228	ierr = NF_INQ_VARID(nid,"psol", varid)
1229	IF (ierr.NE.NF_NOERR) THEN
1230	PRINT*, "Probleme pour avoir identificateur de psol"
1231	CALL abort
1232	ENDIF
1233	#ifdef NC_DOUBLE
1234	ierr = NF_GET_VARA_DOUBLE(nid,varid,debut3,epais3,zx_tmp_2d)
1235	#else
1236	ierr = NF_GET_VARA_REAL(nid,varid,debut3,epais3,zx_tmp_2d)
1237	#endif
1238
1239	IF (ierr.NE.NF_NOERR) THEN
1240	PRINT*, "Probleme pour extraire psol"
1241	CALL abort
1242	ENDIF
1243	CALL gr_fi_lu(1,ngridmx,iim,jjmp1,psol,zx_tmp_2d)
1244	c ----
1245	IF (debug) PRINT *,'lecture SWdnSFC'
1246	ierr = NF_INQ_VARID(nid,"SWdnSFC", varid)
1247	IF (ierr.NE.NF_NOERR) THEN
1248	PRINT*, "Probleme pour avoir identificateur de SWdnSFC"
1249	CALL abort
1250	ENDIF
1251	#ifdef NC_DOUBLE
1252	ierr = NF_GET_VARA_DOUBLE(nid,varid,debut3,epais3,zx_tmp_2d)
1253	#else
1254	ierr = NF_GET_VARA_REAL(nid,varid,debut3,epais3,zx_tmp_2d)
1255	#endif
1256
1257	IF (ierr.NE.NF_NOERR) THEN
1258	PRINT*, "Probleme pour extraire SWdnSFC"
1259	CALL abort
1260	ENDIF
1261	CALL gr_fi_lu(1,ngridmx,iim,jjmp1,solsdn,zx_tmp_2d)
1262	c ----
1263	IF (debug) PRINT *,'lecture SWupSFC'
1264	ierr = NF_INQ_VARID(nid,"SWupSFC", varid)
1265	IF (ierr.NE.NF_NOERR) THEN
1266	PRINT*, "Probleme pour avoir identificateur de SWupSFC"
1267	CALL abort
1268	ENDIF
1269	#ifdef NC_DOUBLE
1270	ierr = NF_GET_VARA_DOUBLE(nid,varid,debut3,epais3,zx_tmp_2d)
1271	#else
1272	ierr = NF_GET_VARA_REAL(nid,varid,debut3,epais3,zx_tmp_2d)
1273	#endif
1274	IF (ierr.NE.NF_NOERR) THEN
1275	PRINT*, "Probleme pour extraire SWupSFC"
1276	CALL abort
1277	ENDIF
1278	CALL gr_fi_lu(1,ngridmx,iim,jjmp1,solsup,zx_tmp_2d)
1279	C
1280	C calcul de l'aledo
1281	DO n=1,ngridmx
1282	IF (solsdn(n).LT.1.0e-3) THEN
1283	albs(n)=0.8
1284	ELSE
1285	albs(n)=solsup(n)/solsdn(n)
1286	ENDIF
1287	ENDDO
1288	c ----
1289	debut4(1) = 1
1290	debut4(2) = 1
1291	debut4(3) = 1
1292	debut4(4) = m
1293	epais4(1) = iim
1294	epais4(2) = jjmp1
1295	epais4(3) = klev
1296	epais4(4) = 1
1297	c
1298	IF (debug) PRINT *,'lecture pres'
1299	ierr = NF_INQ_VARID(nid,"pres", varid)
1300	IF (ierr.NE.NF_NOERR) THEN
1301	PRINT*, "Probleme pour avoir identificateur de pres"
1302	CALL abort
1303	ENDIF
1304	#ifdef NC_DOUBLE
1305	ierr = NF_GET_VARA_DOUBLE(nid,varid,debut4,epais4,zx_tmp_3d)
1306	#else
1307	ierr = NF_GET_VARA_REAL(nid,varid,debut4,epais4,zx_tmp_3d)
1308	#endif
1309	IF (ierr.NE.NF_NOERR) THEN
1310	PRINT*, "Probleme pour extraire pres"
1311	CALL abort
1312	ENDIF
1313	CALL gr_fi_lu(klev,ngridmx,iim,jjmp1,pplay,zx_tmp_3d)
1314	c
1315	DO n = 1, ngridmx
1316	paprs(n,1) = psol(n)
1317	paprs(n,klev+1) = 0.0
1318	ENDDO
1319	DO k = 2, klev
1320	DO n = 1, ngridmx
1321	paprs(n,k) = (pplay(n,k)+pplay(n,k-1))*0.5
1322	ENDDO
1323	ENDDO
1324	c ----
1325	IF (debug) PRINT *,'lecture temp'
1326	ierr = NF_INQ_VARID(nid,"temp", varid)
1327	IF (ierr.NE.NF_NOERR) THEN
1328	PRINT*, "Probleme pour avoir identificateur de temp"
1329	CALL abort
1330	ENDIF
1331	#ifdef NC_DOUBLE
1332	ierr = NF_GET_VARA_DOUBLE(nid,varid,debut4,epais4,zx_tmp_3d)
1333	#else
1334	ierr = NF_GET_VARA_REAL(nid,varid,debut4,epais4,zx_tmp_3d)
1335	#endif
1336	IF (ierr.NE.NF_NOERR) THEN
1337	PRINT*, "Probleme pour extraire temp"
1338	CALL abort
1339	ENDIF
1340	CALL gr_fi_lu(klev,ngridmx,iim,jjmp1,t,zx_tmp_3d)
1341	c ----
1342	IF (debug) PRINT *,'lecture ovap'
1343	ierr = NF_INQ_VARID(nid,"ovap", varid)
1344	IF (ierr.NE.NF_NOERR) THEN
1345	PRINT*, "Probleme pour avoir identificateur de ovap"
1346	CALL abort
1347	ENDIF
1348	#ifdef NC_DOUBLE
1349	ierr = NF_GET_VARA_DOUBLE(nid,varid,debut4,epais4,zx_tmp_3d)
1350	#else
1351	ierr = NF_GET_VARA_REAL(nid,varid,debut4,epais4,zx_tmp_3d)
1352	#endif
1353	IF (ierr.NE.NF_NOERR) THEN
1354	PRINT*, "Probleme pour extraire ovap"
1355	CALL abort
1356	ENDIF
1357	CALL gr_fi_lu(klev,ngridmx,iim,jjmp1,q,zx_tmp_3d)
1358
1359	cozone ----
1360	IF (debug) PRINT *,'ozone'
1361	! ozone(:,:)=0.
1362	ierr = NF_INQ_VARID(nid,"ozone", varid)
1363	IF (ierr.NE.NF_NOERR) THEN
1364	PRINT*, "Probleme pour avoir identificateur d ozone"
1365	CALL abort
1366	ENDIF
1367	! ELSE
1368	#ifdef NC_DOUBLE
1369	ierr = NF_GET_VARA_DOUBLE(nid,varid,debut4,epais4,zx_tmp_3d)
1370	#else
1371	ierr = NF_GET_VARA_REAL(nid,varid,debut4,epais4,zx_tmp_3d)
1372	#endif
1373	IF (ierr.NE.NF_NOERR) THEN
1374	PRINT*, "Probleme pour extraire ozone"
1375	CALL abort
1376	! ENDIF
1377	ENDIF
1378	CALL gr_fi_lu(klev,ngridmx,iim,jjmp1,ozone,zx_tmp_3d)
1379
1380	cozone_daylight
1381	if (reado.eq.2) then
1382	IF (debug) PRINT *,'ozone ozone_daylight'
1383	! ozone(:,:)=0.
1384	ierr = NF_INQ_VARID(nid,"ozone_daylight", varid)
1385	IF (ierr.NE.NF_NOERR) THEN
1386	PRINT*, "Probleme pour avoir identificateur d ozone_dl"
1387	CALL abort
1388	ENDIF
1389	! ELSE
1390	#ifdef NC_DOUBLE
1391	ierr = NF_GET_VARA_DOUBLE(nid,varid,debut4,epais4,zx_tmp_3d)
1392	#else
1393	ierr = NF_GET_VARA_REAL(nid,varid,debut4,epais4,zx_tmp_3d)
1394	#endif
1395	IF (ierr.NE.NF_NOERR) THEN
1396	PRINT*, "Probleme pour extraire ozone_dl"
1397	CALL abort
1398	! ENDIF
1399	ENDIF
1400	CALL gr_fi_lu(klev,ngridmx,iim,jjmp1,ozone_dl,zx_tmp_3d)
1401	endif
1402
1403	c ----
1404	IF (debug) PRINT *,'lecture rneb'
1405	ierr = NF_INQ_VARID(nid,"rneb", varid)
1406	IF (ierr.NE.NF_NOERR) THEN
1407	PRINT*, "Probleme pour avoir identificateur de rneb"
1408	CALL abort
1409	ENDIF
1410	#ifdef NC_DOUBLE
1411	ierr = NF_GET_VARA_DOUBLE(nid,varid,debut4,epais4,zx_tmp_3d)
1412	#else
1413	ierr = NF_GET_VARA_REAL(nid,varid,debut4,epais4,zx_tmp_3d)
1414	#endif
1415	IF (ierr.NE.NF_NOERR) THEN
1416	PRINT*, "Probleme pour extraire rneb"
1417	CALL abort
1418	ENDIF
1419	CALL gr_fi_lu(klev,ngridmx,iim,jjmp1,cldfra,zx_tmp_3d)
1420	c ----
1421	IF (debug) PRINT *,'lecture lwcon'
1422	ierr = NF_INQ_VARID(nid,"lwcon", varid)
1423	IF (ierr.NE.NF_NOERR) THEN
1424	PRINT*, "Probleme pour avoir identificateur de lwcon"
1425	CALL abort
1426	ENDIF
1427	#ifdef NC_DOUBLE
1428	ierr = NF_GET_VARA_DOUBLE(nid,varid,debut4,epais4,zx_tmp_3d)
1429	#else
1430	ierr = NF_GET_VARA_REAL(nid,varid,debut4,epais4,zx_tmp_3d)
1431	#endif
1432	IF (ierr.NE.NF_NOERR) THEN
1433	PRINT*, "Probleme pour extraire lwcon"
1434	CALL abort
1435	ENDIF
1436	CALL gr_fi_lu(klev,ngridmx,iim,jjmp1,cldliq,zx_tmp_3d)
1437	c ----
1438	c Ice water content
1439	IF (debug) PRINT *,'lecture iwcon'
1440	ierr = NF_INQ_VARID(nid,"iwcon", varid)
1441	IF (ierr.NE.NF_NOERR) THEN
1442	PRINT*, "Probleme pour avoir identificateur de iwcon"
1443	CALL abort
1444	ENDIF
1445	#ifdef NC_DOUBLE
1446	ierr = NF_GET_VARA_DOUBLE(nid,varid,debut4,epais4,zx_tmp_3d)
1447	#else
1448	ierr = NF_GET_VARA_REAL(nid,varid,debut4,epais4,zx_tmp_3d)
1449	#endif
1450	IF (ierr.NE.NF_NOERR) THEN
1451	PRINT*, "Probleme pour extraire iwcon"
1452	CALL abort
1453	ENDIF
1454	CALL gr_fi_lu(klev,ngridmx,iim,jjmp1,cldice,zx_tmp_3d)
1455	C
1456	cldwcon=cldliq+cldice
1457	c
1458	ierr = NF_CLOSE(nid)
1459	c
1460	END
1461	C
1462	C ===================================================
1463	C
1464	SUBROUTINE gr_fi_lu(nfield,nlon,iim,jjmp1,fi,lu)
1465	IMPLICIT none
1466	c
1467	c Tranformer une variable de la grille de lecture (hist) a
1468	c la grille physique
1469	C Routine symetrique de gr_fi_ecrit dans physiq
1470	c
1471	INTEGER nfield,nlon,iim,jjmp1, jjm
1472	REAL fi(nlon,nfield), lu(iim*jjmp1,nfield)
1473	c
1474	INTEGER i, n, ig
1475	c
1476	jjm = jjmp1 - 1
1477	DO n = 1, nfield
1478	fi(1,n) = lu(1,n)
1479	fi(nlon,n) = lu(1+jjm*iim,n)
1480	DO ig = 1, nlon - 2
1481	fi(1+ig,n)= lu(iim+ig,n)
1482	ENDDO
1483	ENDDO
1484	RETURN
1485	END
1486	C ---------------------------------------------------------------
1487	SUBROUTINE cstat(ndim,t,ave,std,tmin,tmax,undef)
1488	C
1489	INTEGER ndim, ntot
1490	REAL t(ndim)
1491	REAL ave,std, undef
1492	REAL*8 tsum,tsum2, tave,tstd
1493	REAL tmin,tmax
1494	REAL min, max
1495	C
1496	tsum=0.
1497	ntot=0
1498	tmin=1.E31
1499	tmax=-1.E31
1500	DO n=1,ndim
1501	C print *,t(n),undef
1502	IF (t(n).ne.undef) THEN
1503	tsum=tsum+dble(t(n))
1504	tmin=min(tmin,t(n))
1505	tmax=max(tmax,t(n))
1506	ntot=ntot+1
1507	END IF
1508	END DO
1509	C PRINT *,tsum,ntot
1510	tave=tsum/dble(ntot)
1511	C
1512	tsum2=0.
1513	DO n=1,ndim
1514	IF (t(n).ne.undef) tsum2=tsum2+dble((t(n)-tave)**2)
1515	END DO
1516	tstd=sqrt(tsum2/dble(ntot))
1517	C
1518	ave=sngl(tave)
1519	std=sngl(tstd)
1520	C
1521	END
1522	C
1523	C ---------------------------------------------------------------
1524	SUBROUTINE cminmax(ndim,t,tmin,tmax,undef)
1525	C
1526	IMPLICIT none
1527	INTEGER ndim, ntot, n
1528	REAL t(ndim)
1529	REAL tmin,tmax, undef
1530	REAL min, max
1531	C
1532	tmin=1.E31
1533	tmax=-1.E31
1534	DO n=1,ndim
1535	C print *,t(n),undef
1536	IF (t(n).ne.undef) THEN
1537	tmin=min(tmin,t(n))
1538	tmax=max(tmax,t(n))
1539	END IF
1540	END DO
1541	C
1542	END
1543	C
1544	C ---------------------------------------------------------------
1545	REAL FUNCTION coefcorel(ndim,x,y)
1546	C
1547	INTEGER ndim
1548	REAL x(ndim),y(ndim)
1549	REAL xave, xstd, yave, ystd, covar, xmin, xmax
1550	REAL*8 tsum
1551	C
1552	CALL cstat(ndim, x, xave, xstd, xmin, xmax, undef)
1553	CALL cstat(ndim, y, yave, ystd, xmin, xmax, undef)
1554	C
1555	tsum=0.
1556	DO n=1,ndim
1557	tsum=tsum+dble((x(n)-xave)*(y(n)-yave))
1558	END DO
1559	covar=sngl(tsum/dble(ndim))
1560	coefcorel=covar/(xstd*ystd)
1561	C
1562	END
1563
1564	c-------------------------------------------------------------------
1565	SUBROUTINE cwstat(ndim,t,aire,ave,std,tmin,tmax,undef)
1566	!
1567	INTEGER ndim, ntot
1568	REAL t(ndim),aire(ndim)
1569	REAL ave,std, undef
1570	REAL*8 tsum,tsum2, tave,tstd
1571	REAL*8 asum,aave
1572	REAL tmin,tmax
1573	REAL min, max
1574	!
1575	tsum=0.
1576	asum=0.
1577	ntot=0
1578	tmin=1.E31
1579	tmax=-1.E31
1580	DO n=1,ndim
1581	! print *,t(n),undef
1582	IF (t(n).ne.undef) THEN
1583	tsum=tsum+DBLE(t(n)*aire(n))
1584	asum=asum+DBLE(aire(n))
1585	tmin=min(tmin,t(n))
1586	tmax=max(tmax,t(n))
1587	ntot=ntot+1
1588	END IF
1589	END DO
1590	! PRINT *,tsum,ntot
1591	tave=tsum/asum
1592	aave=asum/dble(ntot)
1593	!
1594	tsum2=0.
1595	DO n=1,ndim
1596	IF (t(n).NE.undef) tsum2=tsum2+DBLE(((t(n)-tave)*2)aire(n)/aave)
1597	END DO
1598	tstd=sqrt(tsum2/dble(ntot))
1599	!
1600	ave=sngl(tave)
1601	std=sngl(tstd)
1602	!
1603	return
1604	!
1605	END
1606

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

Download in other formats:

Original Format