Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

rforcing.F @ 5440

Last change on this file since 5440 was 1547, checked in by idelkadi, 14 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

Rev	Line
[1547]	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