Context Navigation

geo2reg.F90 @ 612

Last change on this file since 612 was 612, checked in by Laurent Fairhead, 20 years ago
Probleme sur les valeurs de w500 qui peuvent etre hors de l'intervalle defini pour la creation des bins LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 18.4 KB

Rev	Line
[606]	1	!
	2	! $Header$
	3	!
	4	!-------------------------------------------------------------------
	5	! Make relationship between any variable and w500 at resolution 2.5x2.5
	6	! -> mean seasonal cycle of the relationship over the period 1985-1989
	7	! -> interannual anomalies (by removing the mean seasonal cycle)
	8	! -> decomposition of interannual anomalies into thermo, dynam and mix components
	9	!
	10	! Specify:
	11	! - undefined value of the data
	12	! - min/max/step of the w500 bins
	13	! - nb of months in timeseries
	14	! - montly index of January 1985
	15	!
	16	! Sandrine Bony, March 2004
	17	! Repris pour faire du netcdf, L. Fairhead 2005/02
	18	!-------------------------------------------------------------------
	19
	20	implicit none
	21
	22	#include "netcdf.inc"
	23
	24	integer imax,jmax,an,mmax,nbinmax,noc
	25	parameter (imax=360,jmax=180,mmax=1700,an=12,nbinmax=200,noc=5)
	26
	27	integer nrun
	28	parameter (nrun=1)
	29
	30	real zmsk(imax,jmax)
	31	real surfmax,zsurf(imax,jmax)
	32	logical ok_msk(imax,jmax)
	33
	34	character*20 namegcm, run
	35
	36	real surftest
	37	real lat0, undef_x
	38	! parameter (lat0 = 30.0) ! limits tropical belt
	39	parameter (undef_x = 1.e+20 ) ! undefined value
	40
	41	real, allocatable, dimension(:,:) :: surf, msk, w, x
	42	real, allocatable, dimension(:,:) :: x_binw
	43	real, allocatable, dimension(:,:) :: w_binw
	44	real, allocatable, dimension(:,:) :: xpdf,nx_binw,sx_binw
	45	real, allocatable, dimension(:) :: xpdfmean
	46
[612]	47	integer nb_bin, nbr_inf, nbr_sup
[606]	48	real min_bin, max_bin, step_bin, w_mult, x_mult
	49	real xpdftot(mmax)
	50	real xpdfmeantot
	51
	52
	53
	54	INTEGER i,j,m,ir,im,jm,lm,itime,n, il, l
	55	real undef, undef_w, pi, x1, w1, msk1
	56	integer :: iostat
	57	character (len=512) :: line_read
	58	character (len=20) :: first_part
	59	character (len=128) :: second_part
	60	character cnat4, bintype4
	61
	62	! netcdf stuff
	63	real :: time
	64	integer, dimension(:), allocatable :: var_dim
	65	INTEGER :: ierr, nvdims, lunout, in_file_id, ierrs, ndims, natts
	66	integer :: lonid, latid, in_time_id, varid, out_file_id, ngatts, levid
	67	integer :: nbin_id, var_bin_id, time_id, var_time_id, in_var_time_id
	68	INTEGER :: var_pdf_id, var_nx_id, var_w_id, var_x_id, outlev_id
	69	integer :: var_lev_id
	70	integer, dimension(4) :: vdims, start, count
	71	character (len=80) :: long_name, varname, attname
	72	character (len=132) :: in_file
	73	real, dimension(:), allocatable ::lat
	74
	75	logical :: var_3d = .false.
	76
	77	integer :: iargc
	78	external iargc
	79	!---------------------------------------------------------------------
	80	!-- General specifications :
	81	!---------------------------------------------------------------------
	82
	83	!***************************** FIN INTERFACE *******************
	84	! Lecture du fichier de configuration
	85	open (20, IOSTAT=iostat, file='config.def',form='formatted')
	86	if (iostat /= 0) THEN
	87	write(,)'Erreur ouverture du fichier config.def'
	88	stop
	89	endif
	90
[611]	91	config: do
	92	read(20,'(A)',iostat=iostat)line_read
	93	if (iostat /= 0) exit
[606]	94	line_read = trim(line_read)
	95	IF (INDEX(line_read, '#') /= 1) THEN
	96	first_part = trim(line_read(1:INDEX(line_read, '=')-1))
	97	second_part = trim(line_read(INDEX(line_read, '=')+1:))
	98	selectcase(first_part)
	99	case('bintype')
	100	bintype = trim(second_part)
	101	case('cnat')
	102	cnat = trim(second_part)
	103	case('lunout')
	104	read(second_part,'(i)') lunout
	105	case('lat0')
	106	read(second_part,'(f)') lat0
	107	end select
	108	endif
[611]	109	enddo config
[606]	110	if (iostat > 0) then
	111	write(lunout,*) &
	112	& 'Probleme de lecture du fichier config.def, iostat = ',iostat
	113	stop
	114	endif
	115	close(20)
	116
	117	! undefined value for the output:
	118	undef = 999.999 ! for output only
	119
	120	!------------------------------------------------------------
	121	! Definition of the w500 bins:
	122	!------------------------------------------------------------
	123
	124	if (bintype.eq.'w500') then
	125	min_bin = -200.0
	126	max_bin = 200.0
	127	step_bin = 10.0
	128	step_bin = 5.0
	129	w_mult = 864.0 ! Pa/s -> hPa/day
	130	! w_mult = 1.0 ! Pa/s -> hPa/day
	131	endif
	132
	133	! scale factor for surface temperature, precip and variable x:
	134	x_mult = 1.0
	135
	136	undef_w = undef
	137
	138	!---------------------------------------------------------------------
	139	!-- preliminaries:
	140	!---------------------------------------------------------------------
	141
	142	if (cnat.ne.'ocea'.and.cnat.ne.'land'.and.cnat.ne.'glob') then
	143	write(,) 'erreur cnat ',cnat
	144	stop
	145	endif
	146
	147	nb_bin = (max_bin-min_bin)/step_bin + 1
	148	if (nb_bin.gt.nbinmax) then
	149	write(,) 'augmenter le nb de bins'
	150	write(,) 'nbinmax, nb_bin = ',nbinmax,nb_bin
	151	stop
	152	endif
	153
	154	! Il faut deux arguments a l appel: fichier et variable
	155	CALL getarg(1, in_file)
	156	if (iargc() == 0 .OR. iargc() /= 2 ) THEN
	157	write(lunout,*)' '
	158	write(lunout,*)' Utilisation de ce programme: '
	159	write(lunout,*)' ./geo2reg nom_de_fichier [variable]'
	160	write(lunout,*) &
	161	& ' ou nom_de_fichier est le nom du fichier a traiter'
	162	write(lunout,*) &
	163	& ' et variable la variable a traiter [optionel]'
	164	write(lunout,*)' '
	165	write(lunout,*)' ./geo2reg -h sort ce message'
	166	write(lunout,*)' '
	167	stop
	168	endif
	169	CALL getarg(2, varname)
	170
	171	! Ouverture du fichier a traiter (histmth)
	172	ierr = NF_OPEN(in_file, NF_NOwrite, in_file_id)
	173	if (ierr /= NF_NOERR) then
	174	write(lunout,*)NF_STRERROR(ierr)
	175	stop
	176	endif
	177	!
	178	! lire im, jm, lm et temps
	179	ierrs = 0
	180	ierr = NF_INQ_DIMID(in_file_id, 'lon', lonid)
	181	ierrs = ierrs + ierr
	182	ierr = NF_INQ_DIMLEN(in_file_id, lonid, im)
	183	ierrs = ierrs + ierr
	184	if (ierrs /= 2 * NF_NOERR) THEN
	185	write(lunout,*)'Pb. avec la lecture de la longitude'
	186	stop
	187	endif
	188	ierrs = 0
	189	ierr = NF_INQ_DIMID(in_file_id, 'lat', latid)
	190	ierrs = ierrs + ierr
	191	ierr = NF_INQ_DIMLEN(in_file_id, latid, jm)
	192	ierrs = ierrs + ierr
	193	allocate (lat(jm))
	194	ierr = NF_INQ_VARID(in_file_id,'lat', latid)
	195	ierrs = ierrs + ierr
	196	ierr = NF_GET_VAR_REAL(in_file_id, latid, lat)
	197	ierrs = ierrs + ierr
	198	if (ierrs /= 4 * NF_NOERR) THEN
	199	write(lunout,*)'Pb. avec la lecture de la latitude'
	200	stop
	201	endif
	202	ierrs = 0
	203	ierr = NF_INQ_DIMID(in_file_id, 'presnivs', levid)
	204	ierrs = ierrs + ierr
	205	ierr = NF_INQ_DIMLEN(in_file_id, levid, lm)
	206	ierrs = ierrs + ierr
	207	if (ierrs /= 2 * NF_NOERR) THEN
	208	write(lunout,*)'Pb. avec la lecture des niveaux verticaux'
	209	stop
	210	endif
	211
	212
	213	ierrs = 0
	214	ierr = NF_INQ_DIMID(in_file_id, 'time_counter', in_time_id)
	215	ierrs = ierrs + ierr
	216	ierr = NF_INQ_DIMLEN(in_file_id, in_time_id, itime)
	217	ierrs = ierrs + ierr
	218	if (ierrs /= 2 * NF_NOERR) THEN
	219	write(lunout,*)'Pb. avec la lecture du temps'
	220	stop
	221	endif
	222
	223	! lecture de l aire des mailles
	224	allocate (surf(im,jm))
	225	ierrs = 0
	226	ierr = NF_INQ_VARID(in_file_id, 'aire', varid)
	227	ierrs = ierrs + ierr
	228	ierr = NF_GET_VAR_REAL(in_file_id, varid, surf)
	229	ierrs = ierrs + ierr
	230	if (ierrs /= 2 * NF_NOERR) THEN
	231	write(lunout,*)'Pb. avec la lecture de l''aire'
	232	stop
	233	endif
	234
	235	!-- compute gridbox areas:
	236
	237
	238	pi = acos(-1.)
	239	surfmax=0.
	240	surftest=0.0
	241	do i = 1, im
	242	do j = 1, jm
	243	surftest=surftest+surf(i,j)
	244	if(surf(i,j).gt.surfmax) surfmax=surf(i,j)
	245	enddo
	246	enddo
[611]	247	write(,) 'sfce totale = 4pi ',surftest/pi
[606]	248
	249	!---------------------------------------------------------------
	250	! Lecture du masque
	251	!
	252	allocate (msk(im,jm))
	253	if (cnat.ne.'glob') then
	254	ierrs = 0
	255	ierr = NF_INQ_VARID(in_file_id, 'pourc_ter', varid)
	256	ierrs = ierrs + ierr
	257	ierr = NF_GET_VAR_REAL(in_file_id, varid, msk)
	258	ierrs = ierrs + ierr
	259	if (ierrs /= 2 * NF_NOERR) THEN
	260	write(lunout,*)'Pb. avec la lecture du masque'
	261	stop
	262	endif
	263	else
	264	do i = 1, im
	265	do j = 1, jm
	266	msk(i,j) = undef
	267	enddo
	268	enddo
	269	endif
	270
	271
	272	do i = 1, im
	273	do j = 1, jm
	274	ok_msk(i,j) = .FALSE.
	275	if (cnat.eq.'ocea' .and. msk(i,j).le.0.5) ok_msk(i,j) = .TRUE.
	276	if (cnat.eq.'land' .and. msk(i,j).gt.0.5) ok_msk(i,j) = .TRUE.
	277	if (cnat.eq.'glob') ok_msk(i,j) = .TRUE.
	278	enddo
	279	enddo
	280
	281	!
	282	! Pour savoir si la variable a traiter est 3D
	283	!
	284	ierrs = 0
	285	ierr = NF_INQ_VARID(in_file_id, varname, varid)
	286	ierrs = ierrs + ierr
	287	ierrs = NF_INQ_VARNDIMS(in_file_id, varid, ndims)
	288	if (ierrs /= 2 * NF_NOERR) THEN
	289	write(lunout,*)'Pb. avec la lecture de ', varname
	290	stop
	291	endif
	292	if (ndims == 4) var_3d = .true.
	293
	294	!--------------------------------------------------------------------
	295	! creation du fichier de sortie et entete
	296	!--------------------------------------------------------------------
	297
	298
	299	ierr = NF_CREATE('out.nc', NF_CLOBBER, out_file_id)
	300	if (ierr /= NF_NOERR) then
	301	write(lunout,*)NF_STRERROR(ierr)
	302	stop
	303	endif
	304	ierr = NF_INQ_NATTS(in_file_id, ngatts)
	305	do i = 1, ngatts
	306	ierr = NF_INQ_ATTNAME(in_file_id, NF_GLOBAL, i, attname)
	307	ierr = NF_COPY_ATT(in_file_id, NF_GLOBAL, attname, out_file_id, NF_GLOBAL)
	308	enddo
	309	!
	310	! definition des dimensions (nbin, level, temps)
	311	! nbin:
	312	ierrs = 0
	313	ierr = NF_DEF_DIM(out_file_id, 'nbin', nb_bin, nbin_id)
	314	ierrs = ierrs + ierr
	315	nvdims = 1
	316	vdims(1) = nbin_id
	317	ierr = NF_DEF_VAR(out_file_id, 'nbin', NF_INT, nvdims, vdims, var_bin_id)
	318	ierrs = ierrs + ierr
	319	ierr = NF_PUT_ATT_TEXT(out_file_id, var_bin_id, 'units', 1, '-')
	320	ierrs = ierrs + ierr
	321	ierr = NF_PUT_ATT_TEXT(out_file_id, var_bin_id, 'long_name', 15, &
	322	& 'number of bins')
	323	ierrs = ierrs + ierr
	324	if (ierrs /= 4 * NF_NOERR) THEN
	325	write(lunout,*)'Pb. dans la definition de nbin'
	326	stop
	327	endif
	328	!
	329	! presnivs
	330	ierrs = 0
	331	ierr = NF_DEF_DIM(out_file_id, 'presnivs', lm, outlev_id)
	332	ierrs = ierrs + ierr
	333	nvdims = 1
	334	vdims(1) = outlev_id
	335	ierr = NF_DEF_VAR(out_file_id,'presnivs',NF_FLOAT, nvdims, vdims, var_lev_id)
	336	ierrs = ierrs + ierr
	337	ierr = NF_INQ_VARID(in_file_id, 'presnivs', levid)
	338	ierrs = ierrs + ierr
	339	ierr = NF_INQ_VARNATTS(in_file_id, levid, natts)
	340	ierrs = ierrs + ierr
	341	do i = 1, natts
	342	attname = ''
	343	ierr = NF_INQ_ATTNAME(in_file_id, levid, i, attname)
	344	ierrs = ierrs + ierr
	345	ierr = NF_COPY_ATT(in_file_id, levid, attname, out_file_id, var_lev_id)
	346	ierrs = ierrs + ierr
	347	enddo
	348	IF (ierrs /= (4 + natts) * NF_NOERR) THEN
	349	write(lunout,*)'Pb. dans la definition de presnivs'
	350	stop
	351	endif
	352
	353	!
	354	! temps:
	355	ierrs = 0
	356	ierr = NF_DEF_DIM(out_file_id, 'time_counter', NF_UNLIMITED, time_id)
	357	ierrs = ierrs + ierr
	358	nvdims = 1
	359	vdims(1) = time_id
	360	ierr = NF_DEF_VAR(out_file_id, 'time_counter', NF_FLOAT, nvdims, vdims, &
	361	& var_time_id)
	362	ierrs = ierrs + ierr
	363	ierr = NF_INQ_VARID(in_file_id, 'time_counter', in_var_time_id)
	364	if (ierr /= NF_NOERR) &
	365	& ierr = NF_INQ_VARID(in_file_id, 't_ave_02592000', in_var_time_id)
	366	if (ierr == NF_NOERR) then
	367	ierr = NF_COPY_ATT(in_file_id, in_var_time_id, 'units', out_file_id, &
	368	& var_time_id)
	369	else
	370	ierr = NF_PUT_ATT_TEXT(out_file_id, var_time_id, 'units', 34, &
	371	& 'seconds since 1860-01-01 00:00:00')
	372	in_var_time_id = 0
	373	endif
	374	ierrs = ierrs + ierr
	375	ierr = NF_PUT_ATT_TEXT(out_file_id, var_time_id, 'long_name', 9, &
	376	& 'Time axis')
	377	ierrs = ierrs + ierr
	378	if (ierrs /= 4 * NF_NOERR) THEN
	379	write(lunout,*)'Pb. dans la definition de time_counter'
	380	stop
	381	endif
	382	!
	383	! Definition des variables a ecrire:
	384	! pdf:
	385	ierrs = 0
	386	nvdims = 2
	387	vdims(1) = nbin_id
	388	vdims(2) = time_id
	389	ierr = NF_DEF_VAR(out_file_id, 'pdf', NF_FLOAT, nvdims, vdims, &
	390	& var_pdf_id)
	391	ierrs = ierrs + ierr
	392	ierr = NF_PUT_ATT_TEXT(out_file_id, var_pdf_id, 'long_name', 16, &
	393	& 'monthly PDF w500')
	394	ierrs = ierrs + ierr
	395	ierr = NF_PUT_ATT_REAL(out_file_id, var_pdf_id, 'missing_value', NF_FLOAT, &
	396	& 1, undef)
	397	ierrs = ierrs + ierr
	398	if (ierrs /= 3 * NF_NOERR) THEN
	399	write(lunout,*)'Pb. dans la definition de pdf'
	400	stop
	401	endif
	402	! nx:
	403	ierrs = 0
	404	nvdims = 2
	405	vdims(1) = nbin_id
	406	vdims(2) = time_id
	407	ierr = NF_DEF_VAR(out_file_id, 'nx', NF_FLOAT, nvdims, vdims, &
	408	& var_nx_id)
	409	ierrs = ierrs + ierr
	410	ierr = NF_PUT_ATT_TEXT(out_file_id, var_nx_id, 'long_name', 16, &
	411	& 'nb of points in w500 bin')
	412	ierrs = ierrs + ierr
	413	if (ierrs /= 2 * NF_NOERR) THEN
	414	write(lunout,*)'Pb. dans la definition de nx'
	415	stop
	416	endif
	417	!
	418	! w:
	419	ierrs = 0
	420	nvdims = 2
	421	vdims(1) = nbin_id
	422	vdims(2) = time_id
	423	ierr = NF_DEF_VAR(out_file_id, 'w', NF_FLOAT, nvdims, vdims, &
	424	& var_w_id)
	425	ierrs = ierrs + ierr
	426	long_name = 'mean relationship w(binw)'
	427	ierr = NF_PUT_ATT_TEXT(out_file_id, var_w_id, 'long_name', 25, &
	428	& long_name)
	429	ierrs = ierrs + ierr
	430	ierr = NF_PUT_ATT_REAL(out_file_id, var_w_id, 'missing_value', NF_FLOAT, &
	431	& 1, undef)
	432	ierrs = ierrs + ierr
	433	if (ierrs /= 3 * NF_NOERR) THEN
	434	write(lunout,*)'Pb. dans la definition de w'
	435	stop
	436	endif
	437	!
	438	! x:
	439	ierrs = 0
	440	nvdims = 2
	441	vdims(1) = nbin_id
	442	vdims(2) = time_id
	443	if (var_3d) then
	444	nvdims = 3
	445	vdims(2) = outlev_id
	446	vdims(3) = time_id
	447	endif
	448	ierr = NF_DEF_VAR(out_file_id, varname, NF_FLOAT, nvdims, vdims, &
	449	& var_x_id)
	450	ierrs = ierrs + ierr
	451	long_name = 'mean relationship '//trim(varname)
	452	ierr = NF_PUT_ATT_TEXT(out_file_id, var_x_id, 'long_name', len(long_name), &
	453	& long_name)
	454	ierrs = ierrs + ierr
	455	ierr = NF_PUT_ATT_REAL(out_file_id, var_x_id, 'missing_value', NF_FLOAT, &
	456	& 1, undef)
	457	ierrs = ierrs + ierr
	458	if (ierrs /= 3 * NF_NOERR) THEN
	459	write(lunout,*)'Pb. dans la definition de ',varname
	460	stop
	461	endif
	462	ierr = NF_ENDDEF(out_file_id)
	463
	464	!
	465	! On passe aux calculs
	466	!
	467	! initialisations:
	468
	469	allocate(xpdfmean(nb_bin))
	470	allocate(xpdf(nb_bin, itime))
	471	allocate(nx_binw(nb_bin, itime))
	472	allocate(sx_binw(nb_bin, itime))
	473	allocate(w_binw(nb_bin, itime))
	474	allocate(x_binw(nb_bin, itime))
	475
	476	!-- temporal loop:
	477	il = 1
	478	if (var_3d) then
	479	il = lm
	480	endif
	481
	482	do l = 1, il
	483	xpdfmeantot = 0.
	484	xpdfmean = 0.0
	485	x_binw = 0.0
	486	xpdf = 0.0
	487	nx_binw = 0.0
	488	sx_binw = 0.0
	489	w_binw = 0.0
	490	xpdftot = 0.
	491	do m = 1, itime ! loop over timeperiod of x_data
	492
	493	! read data:
	494	start(1) = 1
	495	start(2) = 1
	496	start(3) = m
	497	count(1) = im
	498	count(2) = jm
	499	count(3) = 1
	500	! vit. verticale à 500hP
	501	if (.not. allocated(w)) allocate (w(im,jm))
	502	ierrs = 0
	503	ierr = NF_INQ_VARID(in_file_id, 'w500', varid)
	504	ierrs = ierrs + ierr
	505	ierr = NF_GET_VARA_REAL(in_file_id, varid, start, count, w)
	506	ierrs = ierrs + ierr
	507	if (ierrs /= 2* NF_NOERR) THEN
	508	write(lunout,*)'Pb. avec la lecture de la vit. vert. à 500hP'
	509	stop
	510	endif
	511	!
	512	! variable a traiter
	513	if (.not. allocated(x)) allocate (x(im,jm))
	514	if (var_3d) then
	515	start(3) = l
	516	count(3) = 1
	517	start(4) = m
	518	count(4) = 1
	519	endif
	520	ierrs = 0
	521	ierr = NF_INQ_VARID(in_file_id, varname, varid)
	522	ierrs = ierrs + ierr
	523	ierr = NF_GET_VARA_REAL(in_file_id, varid, start, count, x)
	524	ierrs = ierrs + ierr
	525	if (ierrs /= 2 * NF_NOERR) THEN
	526	write(lunout,*)'Pb. avec la lecture de ', varname
	527	stop
	528	endif
	529	!
	530	do j=1,jm
	531	do i=1,im
	532	zmsk(i,j)=0.
	533	enddo
	534	enddo
[612]	535	nbr_inf = 0; nbr_sup = 0
[606]	536	! tropical belt:
	537	do j = 1, jm
	538	if (ABS(lat(j)).le.lat0) then
	539	! loop over grid points:
	540	do i = 1, im
	541	if (x(i,j).ne.undef_x .and. w(i,j).ne.undef_w &
	542	& .and. ok_msk(i,j) ) then
	543	zmsk(i,j)=1.
	544	w1 = w(i,j)*w_mult
	545	x1 = x(i,j)*x_mult
	546	! bin w500:
[612]	547	if (w1 < min_bin) then
	548	ir = 1
	549	nbr_inf = nbr_inf + 1
	550	else if (w1 > max_bin) then
	551	ir = nb_bin
	552	nbr_sup = nbr_sup + 1
	553	else
	554	ir = INT((w1-min_bin)/step_bin) + 1
	555	endif
[606]	556	! monthly PDF :
	557	xpdfmeantot = xpdfmeantot + surf(i,j)
	558	xpdfmean(ir) = xpdfmean(ir) + surf(i,j)
	559	xpdftot(m) = xpdftot(m) + surf(i,j)
	560	xpdf(ir,m) = xpdf(ir,m) + surf(i,j)
	561	! monthly x-w relationship:
	562	nx_binw(ir,m) = nx_binw(ir,m) + 1.0
	563	sx_binw(ir,m) = sx_binw(ir,m) + surf(i,j)
	564	x_binw(ir,m) = x_binw(ir,m) + x1*surf(i,j)
	565	w_binw(ir,m) = w_binw(ir,m) + w1*surf(i,j)
	566
	567	endif
	568	enddo ! i
	569	endif ! lat
	570	enddo ! j
[612]	571	write(lunout,*)'nbre de points ou w500 < ',min_bin,' = ', nbr_inf
	572	write(lunout,*)'nbre de points ou w500 > ',max_bin,' = ', nbr_sup
[606]	573	enddo ! m
	574
	575	! normalize PDF:
	576
	577	do ir = 1, nb_bin
	578	xpdfmean(ir) = xpdfmean(ir)/xpdfmeantot
	579	enddo
	580
	581	do m = 1, itime
	582	if (xpdftot(m).gt.0.) then
	583	do ir = 1, nb_bin
	584	xpdf(ir,m) = xpdf(ir,m)/xpdftot(m)
	585	enddo
	586	else
	587	do ir = 1, nb_bin
	588	xpdf(ir,m) = undef
	589	enddo
	590	endif
	591	enddo
	592
	593
	594	do ir = 1, nb_bin
	595	do m = 1, itime
	596	if (nx_binw(ir,m).gt.1.) then
	597	! if (nx_binw(ir,m).gt.0.) then
	598	x_binw(ir,m) = x_binw(ir,m)/sx_binw(ir,m)
	599	w_binw(ir,m) = w_binw(ir,m)/sx_binw(ir,m)
	600	else
	601	x_binw(ir,m) = undef
	602	w_binw(ir,m) = undef
	603	endif
	604	enddo
	605	enddo ! ir
	606
	607	!
	608	! ecriture du fichier sortie
	609	!
	610
	611	allocate(var_dim(nb_bin))
	612	do ir = 1, nb_bin
	613	var_dim(ir) = min_bin + step_bin * (ir -1)
	614	enddo
	615	ierr = NF_PUT_VAR_INT(out_file_id, var_bin_id, var_dim)
	616	if (ierr /= NF_NOERR) then
	617	write(lunout,*)NF_STRERROR(ierr)
	618	stop
	619	endif
	620	deallocate(var_dim)
	621	!
	622	! ecriture niveaux verticaux
	623	!
	624	allocate(var_dim(lm))
	625	ierrs = 0
	626	ierr = NF_GET_VAR_REAL(in_file_id, levid, var_dim)
	627	ierrs = ierrs + ierr
	628	ierr = NF_PUT_VAR_REAL(out_file_id, var_lev_id, var_dim)
	629	ierrs = ierrs + ierr
	630	if (ierrs /= 2 * NF_NOERR) THEN
	631	write(lunout,*)'Pb. d''ecriture niveaux verticaux'
	632	stop
	633	endif
[612]	634	deallocate(var_dim)
[606]	635
	636
	637
	638	do m = 1, itime
	639	start(1) = 1
	640	start(2) = m
	641	count(1) = nb_bin
	642	count(2) = 1
	643	ierrs = 0
	644	!
	645	! ecriture temps:
	646	if (in_var_time_id == 0) then
	647	time = (m - 1) * 86400. * 30.
	648	else
	649	ierr = NF_GET_VAR1_REAL(in_file_id, in_var_time_id, m, time)
	650	endif
	651	ierr = NF_PUT_VAR1_REAL(out_file_id, var_time_id, m, time)
	652	ierrs = ierrs + ierr
	653	!
	654	! pdf, nb_bin, x, w:
	655	ierr = NF_PUT_VARA_REAL(out_file_id, var_pdf_id, start, count, &
	656	& xpdf(1,m))
	657	ierrs = ierrs + ierr
	658	ierr = NF_PUT_VARA_REAL(out_file_id, var_nx_id, start, count, &
	659	& nx_binw(1,m))
	660	ierrs = ierrs + ierr
	661	ierr = NF_PUT_VARA_REAL(out_file_id, var_w_id, start, count, &
	662	& w_binw(1,m))
	663	ierrs = ierrs + ierr
	664	if (var_3d) then
	665	start(2) = l
	666	count(2) = 1
	667	start(3) = m
	668	count(4) = 1
	669	endif
	670	ierr = NF_PUT_VARA_REAL(out_file_id, var_x_id, start, count, &
	671	& x_binw(1,m))
	672	ierrs = ierrs + ierr
	673	if (ierrs /= 5 * NF_NOERR) THEN
	674	write(lunout,*)'Pb. d''ecriture'
	675	stop
	676	endif
	677	enddo ! m
	678	enddo ! l
	679
	680	ierr = NF_CLOSE(out_file_id)
	681
	682	end

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

Download in other formats:

Original Format