Context Navigation

wstats_mod.F90 @ 3567

Last change on this file since 3567 was 3282, checked in by emillour, 9 months ago
Mars PCM: Fix a buggy behavior of wstats: it would fail when using a stats.def file and when the first (at every time step) variable sent to wstats() was not part of those listed in stats.def EM
File size: 28.9 KB

Line
1	module wstats_mod
2
3	! module containing parameters and routines to generate the "stats.nc" file
4	! which will contain a mean statistical day of variables extracted at set
5	! times of day every day of the simulation, and the standard deviations thereof
6
7	implicit none
8
9	logical,save :: callstats ! global flag to trigger generating a stats.nc
10	! file or not. Initialized in conf_gcm()
11	!$OMP THREADPRIVATE(callstats)
12
13	integer,save :: istats ! calculate stats every istats physics timestep,
14	! starting at first call. Initialized by inistats()
15	!$OMP THREADPRIVATE(istats)
16
17	integer,parameter :: istime=12 ! number of time steps per sol to store
18
19	integer,save :: count(istime) ! count tab to know the variable record
20	!$OMP THREADPRIVATE(count)
21
22	contains
23
24	subroutine wstats(ngrid,nom,titre,unite,dim,px)
25
26	! Main routine to call to trigger writing a given field to the "stats.nc" file
27
28	use mod_phys_lmdz_para, only : is_mpi_root, is_master, klon_mpi_begin
29	use mod_phys_lmdz_transfert_para, only: gather, bcast
30	use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo, &
31	nbp_lon, nbp_lat, nbp_lev, &
32	grid_type, unstructured
33	implicit none
34
35	include "netcdf.inc"
36
37	integer,intent(in) :: ngrid
38	character (len=*),intent(in) :: nom,titre,unite
39	integer,intent(in) :: dim
40	real,intent(in) :: px(ngrid,nbp_lev)
41	real,allocatable,save :: mean3d(:,:,:),sd3d(:,:,:),dx3(:,:,:)
42	real,allocatable,save :: mean2d(:,:),sd2d(:,:),dx2(:,:)
43	character (len=50) :: namebis
44	character (len=50), save :: firstvar
45	!$OMP THREADPRIVATE(mean3d,sd3d,dx3,mean2d,sd2d,dx2,firstvar)
46	integer :: ierr,varid,nbdim,nid
47	integer :: meanid,sdid
48	integer, dimension(4) :: id,start,sizes
49	logical, save :: firstcall=.TRUE.
50	integer,save :: indx
51	integer, save :: step=0
52	!$OMP THREADPRIVATE(firstcall,indx,step)
53	integer :: l,i,j,ig0,n
54
55	! Added to read an optional stats.def file to select outputs
56	logical,save :: check_stats_def=.TRUE. ! to keep track of if stats.def was checked
57	logical,save :: stats_def ! .true. if there is a stats.def file
58	integer,save :: n_name_def ! number of fields to output in stats.nc
59	integer,parameter :: n_name_def_max=199 ! max number of fields to output
60	character(len=120),save :: name_def(n_name_def_max)
61	logical :: getout ! to trigger an early exit if variable not in output list
62
63	!$OMP THREADPRIVATE(check_stats_def,stats_def,n_name_def,name_def)
64
65	! Added to work in parallel mode
66	#ifdef CPP_PARA
67	real px3_glop(klon_glo,nbp_lev) ! to store a 3D data set on global physics grid
68	real px3_glo(nbp_lon,nbp_lat,nbp_lev) ! to store a global 3D data set on global lonxlat grid
69	real px2_glop(klon_glo) ! to store a 2D data set on global physics grid
70	real px2_glo(nbp_lon,nbp_lat) ! to store a 2D data set on global lonxlat grid
71	real px2(ngrid)
72	real px3(ngrid,nbp_lev)
73	#else
74	! When not running in parallel mode:
75	real px3_glop(ngrid,nbp_lev) ! to store a 3D data set on global physics grid
76	real px3_glo(nbp_lon,nbp_lat,nbp_lev) ! to store a global 3D data set on global lonxlat grid
77	real px2_glop(ngrid) ! to store a 2D data set on global physics grid
78	real px2_glo(nbp_lon,nbp_lat) ! to store a 2D data set on global lonxlat grid
79	#endif
80
81	! 0. Preliminary stuff
82	if (callstats.eqv..false.) then
83	! exit because creating/writing stats.nc not requested by user
84	return
85	endif
86
87	if (grid_type==unstructured) then
88	! exit because non-structured grid case is not handled
89	return
90	endif
91
92	if (check_stats_def) then
93	! this step only needs be made a single unique time
94	! hence the chack_stats_def flag.
95
96	if (is_master) then
97	! Only the master reads the file
98	! open stats.def definition file if there is one
99	open(99,file="stats.def",status='old',form='formatted',&
100	iostat=ierr)
101	if (ierr.eq.0) then
102	stats_def=.true. ! yes there is a stats.def file
103	write(,) "*****************"
104	write(,) "Reading stats.def"
105	write(,) "*****************"
106	do n=1,n_name_def_max
107	read(99,fmt='(a)',end=88) name_def(n)
108	write(,) 'Output in stats: ', trim(name_def(n))
109	enddo
110	88 continue
111	! check there is no overflow
112	if (n.ge.n_name_def_max) then
113	write(,) "n_name_def_max too small in wstats:",n
114	call abort_physic("wstats","n_name_def_max too small",1)
115	endif
116	n_name_def=n-1
117	close(99)
118	else
119	stats_def=.false. ! no stats.def file; output all fields sent to wstats
120	endif ! of if (ierr.eq.0)
121
122	endif ! of if (is_master)
123
124	! share the info gathered by master with all other cores
125	call bcast(stats_def)
126	call bcast(n_name_def)
127	do n=1,n_name_def
128	call bcast(name_def(n))
129	enddo
130
131	check_stats_def=.false.
132	endif ! of if (check_stats_def)
133
134	! 1. Initialization (creation of stats.nc file)
135	if (firstcall) then
136	firstcall=.false.
137
138	! identify the first variable name with which wstats is called in
139	! a time step
140	if (stats_def) then
141	! check if "nom" is in the wanted output list. If it is then
142	! it is the sought "first variable"; if not then skip all the
143	! rest and reset "firstcall" to true
144	getout=.true.
145	do n=1,n_name_def
146	! look for the variable's name in the list
147	if (trim(name_def(n)).eq.nom) then
148	getout=.false.
149	! found it, no need to scan the rest of the list exit loop
150	firstvar=trim((nom))
151	exit
152	endif
153	enddo
154	if (getout) then
155	! variable not in the list so exit routine now
156	firstcall=.true.
157	return
158	endif
159	else
160	! without a stats.def, the very first call to stats contains first
161	! variable name
162	firstvar=trim((nom))
163	endif ! of if (stats_def)
164
165	! initialize stats.nc file
166	call inistats(ierr)
167
168	! allocate arrays that will be needed later
169	if (klon_glo>1) then ! general case, 3D GCM
170	allocate(mean3d(nbp_lon+1,nbp_lat,nbp_lev))
171	allocate(sd3d(nbp_lon+1,nbp_lat,nbp_lev))
172	allocate(dx3(nbp_lon+1,nbp_lat,nbp_lev))
173	allocate(mean2d(nbp_lon+1,nbp_lat))
174	allocate(sd2d(nbp_lon+1,nbp_lat))
175	allocate(dx2(nbp_lon+1,nbp_lat))
176	else ! 1D model
177	allocate(mean3d(1,1,nbp_lev))
178	allocate(sd3d(1,1,nbp_lev))
179	allocate(dx3(1,1,nbp_lev))
180	allocate(mean2d(1,1))
181	allocate(sd2d(1,1))
182	allocate(dx2(1,1))
183	endif
184
185	endif ! of if (firstcall)
186
187	if (firstvar==nom) then ! If we're back to the first variable, increment time counter
188	step = step + 1
189	endif
190
191	if (mod(step,istats).ne.0) then
192	! if its not time to write to file, exit
193	RETURN
194	endif
195
196	! Exit if there is a stats.def file and the variable is not in the list
197	if (stats_def) then
198	getout=.true.
199	do n=1,n_name_def
200	! look for the variable's name in the list
201	if (trim(name_def(n)).eq.nom) then
202	getout=.false.
203	! found it, no need to scan the rest of the list exit loop
204	exit
205	endif
206	enddo
207	if (getout) then
208	! variable not in the list so exit routine now
209	return
210	endif
211	endif ! of if (stats_def)
212
213	! collect fields on a global physics grid
214	#ifdef CPP_PARA
215	if (dim.eq.3) then
216	px3(1:ngrid,1:nbp_lev)=px(1:ngrid,1:nbp_lev)
217	! Gather fieds on a "global" (without redundant longitude) array
218	call Gather(px3,px3_glop)
219	!$OMP MASTER
220	if (is_mpi_root) then
221	call Grid1Dto2D_glo(px3_glop,px3_glo)
222	if (klon_glo>1) then ! general case (unless in 1D)
223	! copy dx3_glo() to dx3(:) and add redundant longitude
224	dx3(1:nbp_lon,:,:)=px3_glo(1:nbp_lon,:,:)
225	dx3(nbp_lon+1,:,:)=dx3(1,:,:)
226	endif
227	endif
228	!$OMP END MASTER
229	!$OMP BARRIER
230	else ! dim.eq.2
231	! Gather fieds on a "global" (without redundant longitude) array
232	px2(:)=px(:,1)
233	call Gather(px2,px2_glop)
234	!$OMP MASTER
235	if (is_mpi_root) then
236	call Grid1Dto2D_glo(px2_glop,px2_glo)
237	if (klon_glo>1) then ! general case (unless in 1D)
238	! copy px2_glo() to dx2(:) and add redundant longitude
239	dx2(1:nbp_lon,:)=px2_glo(1:nbp_lon,:)
240	dx2(nbp_lon+1,:)=dx2(1,:)
241	endif
242	endif
243	!$OMP END MASTER
244	!$OMP BARRIER
245	endif
246	#else
247	if (dim.eq.3) then
248	px3_glop(:,1:nbp_lev)=px(:,1:nbp_lev)
249	! Passage variable physique --> variable dynamique
250	DO l=1,nbp_lev
251	DO i=1,nbp_lon
252	px3_glo(i,1,l)=px(1,l)
253	px3_glo(i,nbp_lat,l)=px(ngrid,l)
254	ENDDO
255	DO j=2,nbp_lat-1
256	ig0= 1+(j-2)*nbp_lon
257	DO i=1,nbp_lon
258	px3_glo(i,j,l)=px(ig0+i,l)
259	ENDDO
260	ENDDO
261	ENDDO
262	else ! dim.eq.2
263	px2_glop(:)=px(:,1)
264	! Passage variable physique --> physique dynamique
265	DO i=1,nbp_lon
266	px2_glo(i,1)=px(1,1)
267	px2_glo(i,nbp_lat)=px(ngrid,1)
268	ENDDO
269	DO j=2,nbp_lat-1
270	ig0= 1+(j-2)*nbp_lon
271	DO i=1,nbp_lon
272	px2_glo(i,j)=px(ig0+i,1)
273	ENDDO
274	ENDDO
275	endif
276	#endif
277
278	! 2. Write field to file
279
280	if (is_master) then
281	! only master needs do this
282
283	ierr = NF_OPEN("stats.nc",NF_WRITE,nid)
284
285	namebis=trim(nom)
286
287	! test: check if that variable already exists in file
288	ierr= NF_INQ_VARID(nid,namebis,meanid)
289
290	if (ierr.ne.NF_NOERR) then
291	! variable not in file, create/define it
292	if (firstvar==nom) then
293	indx=1
294	count(:)=0
295	endif
296
297
298	!declaration de la variable
299
300	! choix du nom des coordonnees
301	ierr= NF_INQ_DIMID(nid,"longitude",id(1))
302	ierr= NF_INQ_DIMID(nid,"latitude",id(2))
303	if (dim.eq.3) then
304	ierr= NF_INQ_DIMID(nid,"altitude",id(3))
305	ierr= NF_INQ_DIMID(nid,"Time",id(4))
306	nbdim=4
307	else if (dim==2) then
308	ierr= NF_INQ_DIMID(nid,"Time",id(3))
309	nbdim=3
310	endif
311
312	write (,) "====================="
313	write (,) "STATS: creating ",nom
314	namebis=trim(nom)
315	call def_var_stats(nid,namebis,titre,unite,nbdim,id,meanid,ierr)
316	if (dim.eq.3) then
317	call inivar(nid,meanid,size(px3_glop,1),dim,indx,px3_glop,ierr)
318	else ! dim.eq.2
319	call inivar(nid,meanid,size(px2_glop,1),dim,indx,px2_glop,ierr)
320	endif
321	namebis=trim(nom)//"_sd"
322	call def_var_stats(nid,namebis,trim(titre)//" total standard deviation over the season",unite,nbdim,id,sdid,ierr)
323	if (dim.eq.3) then
324	call inivar(nid,sdid,size(px3_glop,1),dim,indx,px3_glop,ierr)
325	else ! dim.eq.2
326	call inivar(nid,sdid,size(px2_glop,1),dim,indx,px2_glop,ierr)
327	endif
328
329	ierr= NF_CLOSE(nid)
330	return
331
332	else
333	! variable found in file
334	namebis=trim(nom)//"_sd"
335	ierr= NF_INQ_VARID(nid,namebis,sdid)
336
337	endif
338
339	if (firstvar==nom) then
340	count(indx)=count(int(indx))+1
341	indx=indx+1
342	if (indx>istime) then
343	indx=1
344	endif
345	endif
346
347	if (count(indx)==0) then
348	! very first time we write the variable to file
349	if (dim.eq.3) then
350	start=(/1,1,1,indx/)
351	if (klon_glo>1) then !general case
352	sizes=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
353	else
354	sizes=(/1,1,nbp_lev,1/)
355	endif
356	mean3d(:,:,:)=0
357	sd3d(:,:,:)=0
358	else if (dim.eq.2) then
359	start=(/1,1,indx,0/)
360	if (klon_glo>1) then !general case
361	sizes=(/nbp_lon+1,nbp_lat,1,0/)
362	else
363	sizes=(/1,1,1,0/)
364	endif
365	mean2d(:,:)=0
366	sd2d(:,:)=0
367	endif
368	else
369	! load values from file
370	if (dim.eq.3) then
371	start=(/1,1,1,indx/)
372	if (klon_glo>1) then !general case
373	sizes=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
374	else ! 1D model case
375	sizes=(/1,1,nbp_lev,1/)
376	endif
377	#ifdef NC_DOUBLE
378	ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,sizes,mean3d)
379	ierr = NF_GET_VARA_DOUBLE(nid,sdid,start,sizes,sd3d)
380	#else
381	ierr = NF_GET_VARA_REAL(nid,meanid,start,sizes,mean3d)
382	ierr = NF_GET_VARA_REAL(nid,sdid,start,sizes,sd3d)
383	#endif
384	if (ierr.ne.NF_NOERR) then
385	write (,) "wstats error reading :",trim(nom)
386	write (,) NF_STRERROR(ierr)
387	call abort_physic("wstats","Failed reading "//trim(nom),1)
388	endif
389
390	else if (dim.eq.2) then
391	start=(/1,1,indx,0/)
392	if (klon_glo>1) then ! general case
393	sizes=(/nbp_lon+1,nbp_lat,1,0/)
394	else
395	sizes=(/1,1,1,0/)
396	endif
397	#ifdef NC_DOUBLE
398	ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,sizes,mean2d)
399	ierr = NF_GET_VARA_DOUBLE(nid,sdid,start,sizes,sd2d)
400	#else
401	ierr = NF_GET_VARA_REAL(nid,meanid,start,sizes,mean2d)
402	ierr = NF_GET_VARA_REAL(nid,sdid,start,sizes,sd2d)
403	#endif
404	if (ierr.ne.NF_NOERR) then
405	write (,) "wstats error reading :",trim(nom)
406	write (,) NF_STRERROR(ierr)
407	call abort_physic("wstats","Failed reading "//trim(nom),1)
408	endif
409	endif
410	endif ! of if (count(indx)==0)
411
412
413	! 2.1. Build dx* (data on lon-lat grid, with redundant longitude)
414
415	if (dim.eq.3) then
416	dx3(1:nbp_lon,:,:)=px3_glo(:,:,:)
417	IF (klon_glo>1) THEN ! in 3D, add redundant longitude point
418	dx3(nbp_lon+1,:,:)=dx3(1,:,:)
419	ENDIF
420	else ! dim.eq.2
421	dx2(1:nbp_lon,:)=px2_glo(:,:)
422	IF (klon_glo>1) THEN ! in 3D, add redundant longitude point
423	dx2(nbp_lon+1,:)=dx2(1,:)
424	ENDIF
425	endif
426
427
428	! 2.2. Add current values to previously stored sums
429
430	if (dim.eq.3) then
431
432	mean3d(:,:,:)=mean3d(:,:,:)+dx3(:,:,:)
433	sd3d(:,:,:)=sd3d(:,:,:)+dx3(:,:,:)**2
434
435	#ifdef NC_DOUBLE
436	ierr = NF_PUT_VARA_DOUBLE(nid,meanid,start,sizes,mean3d)
437	ierr = NF_PUT_VARA_DOUBLE(nid,sdid,start,sizes,sd3d)
438	#else
439	ierr = NF_PUT_VARA_REAL(nid,meanid,start,sizes,mean3d)
440	ierr = NF_PUT_VARA_REAL(nid,sdid,start,sizes,sd3d)
441	#endif
442	if (ierr.ne.NF_NOERR) then
443	write (,) "wstats error writing :",trim(nom)
444	write (,) NF_STRERROR(ierr)
445	call abort_physic("wstats","Failed writing "//trim(nom),1)
446	endif
447
448	else if (dim.eq.2) then
449
450	mean2d(:,:)= mean2d(:,:)+dx2(:,:)
451	sd2d(:,:)=sd2d(:,:)+dx2(:,:)**2
452
453	#ifdef NC_DOUBLE
454	ierr = NF_PUT_VARA_DOUBLE(nid,meanid,start,sizes,mean2d)
455	ierr = NF_PUT_VARA_DOUBLE(nid,sdid,start,sizes,sd2d)
456	#else
457	ierr = NF_PUT_VARA_REAL(nid,meanid,start,sizes,mean2d)
458	ierr = NF_PUT_VARA_REAL(nid,sdid,start,sizes,sd2d)
459	#endif
460	if (ierr.ne.NF_NOERR) then
461	write (,) "wstats error writing :",trim(nom)
462	write(,) "start:",start
463	write(,) "sizes:",sizes
464	write(,) "mean2d:",mean2d
465	write(,) "sd2d:",sd2d
466	write (,) NF_STRERROR(ierr)
467	call abort_physic("wstats","Failed writing "//trim(nom),1)
468	endif
469
470	endif ! of if (dim.eq.3) elseif (dim.eq.2)
471
472	ierr= NF_CLOSE(nid)
473	endif ! of if (is_master)
474
475	end subroutine wstats
476
477	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
478
479	subroutine inistats(ierr)
480
481	! routine to initialize the stats file (i.e. create file, write
482	! time independent coordinates, etc.)
483
484	use mod_phys_lmdz_para, only : is_master
485	USE vertical_layers_mod, ONLY: ap,bp,aps,bps,preff,pseudoalt,presnivs
486	USE nrtype, ONLY: pi
487	USE time_phylmdz_mod, ONLY: daysec,dtphys
488	USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
489	USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev
490	implicit none
491
492	include "netcdf.inc"
493
494	integer,intent(out) :: ierr
495
496	integer :: nid
497	integer :: l,nsteppd
498	real, dimension(nbp_lev) :: sig_s
499	real,allocatable :: lon_reg_ext(:) ! extended longitudes
500	integer :: idim_lat,idim_lon,idim_llm,idim_llmp1,idim_time
501	real, dimension(istime) :: lt
502	integer :: nvarid
503
504
505	IF (nbp_lon*nbp_lat==1) THEN
506	! 1D model
507	ALLOCATE(lon_reg_ext(1))
508	ELSE
509	! 3D model
510	ALLOCATE(lon_reg_ext(nbp_lon+1))
511	ENDIF
512
513	write (,)
514	write (,) ' \|\| STATS \|\|'
515	write (,)
516	write (,) 'daysec',daysec
517	write (,) 'dtphys',dtphys
518	nsteppd=nint(daysec/dtphys)
519	write (,) 'nsteppd=',nsteppd
520
521	if (abs(float(nsteppd)-daysec/dtphys).gt.1.e-8*daysec) then
522	call abort_physic("inistats",'1 sol .ne. n physics steps',1)
523	endif
524
525	if(mod(nsteppd,istime).ne.0) then
526	call abort_physic("inistats",'1 sol .ne. n*istime physics steps',1)
527	endif
528
529	istats=nsteppd/istime
530	write (,) 'istats=',istats
531	write (,) 'Storing ',istime,'times per day'
532	write (,) 'thus every ',istats,'physical timestep '
533	write (,)
534
535	do l= 1, nbp_lev
536	sig_s(l)=((ap(l)+ap(l+1))/preff+bp(l)+bp(l+1))/2.
537	pseudoalt(l)=-10.*log(presnivs(l)/preff)
538	enddo
539
540	lon_reg_ext(1:nbp_lon)=lon_reg(1:nbp_lon)
541	IF (nbp_lon*nbp_lat/=1) THEN
542	! In 3D, add extra redundant point (180 degrees,
543	! since lon_reg starts at -180)
544	lon_reg_ext(nbp_lon+1)=-lon_reg_ext(1)
545	ENDIF
546
547	if (is_master) then
548	! only the master needs do this
549	ierr = NF_CREATE("stats.nc",IOR(NF_CLOBBER,NF_64BIT_OFFSET),nid)
550	if (ierr.ne.NF_NOERR) then
551	write (,) NF_STRERROR(ierr)
552	call abort_physic("inistats",'failed creating stats.nc',1)
553	endif
554
555	ierr = NF_DEF_DIM (nid, "latitude", nbp_lat, idim_lat)
556	IF (nbp_lon*nbp_lat==1) THEN
557	ierr = NF_DEF_DIM (nid, "longitude", 1, idim_lon)
558	ELSE
559	ierr = NF_DEF_DIM (nid, "longitude", nbp_lon+1, idim_lon)
560	ENDIF
561	ierr = NF_DEF_DIM (nid, "altitude", nbp_lev, idim_llm)
562	ierr = NF_DEF_DIM (nid, "llmp1", nbp_lev+1, idim_llmp1)
563	ierr = NF_DEF_DIM (nid, "Time", NF_UNLIMITED, idim_time)
564
565	ierr = NF_ENDDEF(nid)
566
567	call def_var_stats(nid,"Time","Time", &
568	"days since 0000-00-0 00:00:00",1, &
569	[idim_time],nvarid,ierr)
570	! Time is initialised later by mkstats subroutine
571
572	call def_var_stats(nid,"latitude","latitude", &
573	"degrees_north",1,[idim_lat],nvarid,ierr)
574	#ifdef NC_DOUBLE
575	ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lat_reg/pi*180)
576	#else
577	ierr = NF_PUT_VAR_REAL (nid,nvarid,lat_reg/pi*180)
578	#endif
579
580	call def_var_stats(nid,"longitude","East longitude", &
581	"degrees_east",1,[idim_lon],nvarid,ierr)
582	#ifdef NC_DOUBLE
583	ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lon_reg_ext/pi*180)
584	#else
585	ierr = NF_PUT_VAR_REAL (nid,nvarid,lon_reg_ext/pi*180)
586	#endif
587
588	! Niveaux verticaux, aps et bps
589	ierr = NF_REDEF (nid)
590	#ifdef NC_DOUBLE
591	ierr = NF_DEF_VAR (nid,"altitude", NF_DOUBLE, 1,[idim_llm],nvarid)
592	#else
593	ierr = NF_DEF_VAR (nid,"altitude", NF_FLOAT, 1,[idim_llm],nvarid)
594	#endif
595	ierr = NF_PUT_ATT_TEXT (nid,nvarid,"long_name",8,"altitude")
596	ierr = NF_PUT_ATT_TEXT (nid,nvarid,'units',2,"km")
597	ierr = NF_PUT_ATT_TEXT (nid,nvarid,'positive',2,"up")
598	ierr = NF_ENDDEF(nid)
599	#ifdef NC_DOUBLE
600	ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,pseudoalt)
601	#else
602	ierr = NF_PUT_VAR_REAL (nid,nvarid,pseudoalt)
603	#endif
604	call def_var_stats(nid,"aps","hybrid pressure at midlayers", &
605	" ",1,[idim_llm],nvarid,ierr)
606	#ifdef NC_DOUBLE
607	ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,aps)
608	#else
609	ierr = NF_PUT_VAR_REAL (nid,nvarid,aps)
610	#endif
611
612	call def_var_stats(nid,"bps","hybrid sigma at midlayers", &
613	" ",1,[idim_llm],nvarid,ierr)
614	#ifdef NC_DOUBLE
615	ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,bps)
616	#else
617	ierr = NF_PUT_VAR_REAL (nid,nvarid,bps)
618	#endif
619
620	ierr=NF_CLOSE(nid)
621
622	endif ! of if (is_master)
623
624	end subroutine inistats
625
626	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
627
628	subroutine inivar(nid,varid,ngrid,dim,indx,px,ierr)
629
630	! routine to initialize writing a variable in the stats file
631
632	use mod_grid_phy_lmdz, only : nbp_lon, nbp_lat, nbp_lev, klon_glo
633
634	implicit none
635
636	include "netcdf.inc"
637
638	integer, intent(in) :: nid,varid,dim,indx,ngrid
639	real, dimension(ngrid,nbp_lev), intent(in) :: px
640	integer, intent(out) :: ierr
641
642	integer :: l,i,j,ig0
643	integer, dimension(4) :: start,sizes
644	real, dimension(nbp_lon+1,nbp_lat,nbp_lev) :: dx3
645	real, dimension(nbp_lon+1,nbp_lat) :: dx2
646	real :: dx3_1d(nbp_lev) ! for 1D outputs
647	real :: dx2_1d ! for 1D outputs
648
649	if (dim.eq.3) then
650
651	start=(/1,1,1,indx/)
652	if (klon_glo>1) then ! general 3D case
653	sizes=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
654	else
655	sizes=(/1,1,nbp_lev,1/)
656	endif
657
658	! Passage variable physique --> variable dynamique
659
660	if (klon_glo>1) then ! general case
661	DO l=1,nbp_lev
662	DO i=1,nbp_lon+1
663	dx3(i,1,l)=px(1,l)
664	dx3(i,nbp_lat,l)=px(ngrid,l)
665	ENDDO
666	DO j=2,nbp_lat-1
667	ig0= 1+(j-2)*nbp_lon
668	DO i=1,nbp_lon
669	dx3(i,j,l)=px(ig0+i,l)
670	ENDDO
671	dx3(nbp_lon+1,j,l)=dx3(1,j,l)
672	ENDDO
673	ENDDO
674	else ! 1D model case
675	dx3_1d(1:nbp_lev)=px(1,1:nbp_lev)
676	endif
677
678	#ifdef NC_DOUBLE
679	if (klon_glo>1) then
680	ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,sizes,dx3)
681	else
682	ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,sizes,dx3_1d)
683	endif
684	#else
685	if (klon_glo>1) then
686	ierr = NF_PUT_VARA_REAL(nid,varid,start,sizes,dx3)
687	else
688	ierr = NF_PUT_VARA_REAL(nid,varid,start,sizes,dx3_1d)
689	endif
690	#endif
691	if (ierr.ne.NF_NOERR) then
692	write (,) "inivar error writing variable"
693	write (,) NF_STRERROR(ierr)
694	call abort_physic("inivar","error writing variable",1)
695	endif
696
697	else if (dim.eq.2) then
698
699	start=(/1,1,indx,0/)
700	if (klon_glo>1) then ! general 3D case
701	sizes=(/nbp_lon+1,nbp_lat,1,0/)
702	else
703	sizes=(/1,1,1,0/)
704	endif
705
706	! Passage variable physique --> physique dynamique
707
708	if (klon_glo>1) then ! general case
709	DO i=1,nbp_lon+1
710	dx2(i,1)=px(1,1)
711	dx2(i,nbp_lat)=px(ngrid,1)
712	ENDDO
713	DO j=2,nbp_lat-1
714	ig0= 1+(j-2)*nbp_lon
715	DO i=1,nbp_lon
716	dx2(i,j)=px(ig0+i,1)
717	ENDDO
718	dx2(nbp_lon+1,j)=dx2(1,j)
719	ENDDO
720	else ! 1D model case
721	dx2_1d=px(1,1)
722	endif
723
724	#ifdef NC_DOUBLE
725	if (klon_glo>1) then
726	ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,sizes,dx2)
727	else
728	ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,sizes,dx2_1d)
729	endif
730	#else
731	if (klon_glo>1) then
732	ierr = NF_PUT_VARA_REAL(nid,varid,start,sizes,dx2)
733	else
734	ierr = NF_PUT_VARA_REAL(nid,varid,start,sizes,dx2_1d)
735	endif
736	#endif
737	if (ierr.ne.NF_NOERR) then
738	write (,) "inivar error writing variable"
739	write (,) NF_STRERROR(ierr)
740	call abort_physic("inivar","error writing variable",1)
741	endif
742
743	endif
744
745	end subroutine inivar
746
747	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
748
749	subroutine def_var_stats(nid,name,title,units,nbdim,dimids,nvarid,ierr)
750
751	! This subroutine defines variable 'name' in a (pre-existing and opened)
752	! NetCDF file (known from its NetCDF ID 'nid').
753	! The number of dimensions 'nbdim' of the variable, as well as the IDs of
754	! corresponding dimensions must be set (in array 'dimids').
755	! Upon successfull definition of the variable, 'nvarid' contains the
756	! NetCDF ID of the variable.
757	! The variables' attributes 'title' (Note that 'long_name' would be more
758	! appropriate) and 'units' are also set.
759
760	implicit none
761
762	include "netcdf.inc"
763
764	integer,intent(in) :: nid ! NetCDF file ID
765	character(len=*),intent(in) :: name ! the variable's name
766	character(len=*),intent(in) :: title ! 'title' attribute of variable
767	character(len=*),intent(in) :: units ! 'units' attribute of variable
768	integer,intent(in) :: nbdim ! number of dimensions of the variable
769	integer,dimension(nbdim),intent(in) :: dimids ! NetCDF IDs of the dimensions
770	! the variable is defined along
771	integer,intent(out) :: nvarid ! NetCDF ID of the variable
772	integer,intent(out) :: ierr ! returned NetCDF staus code
773
774	! 1. Switch to NetCDF define mode
775	ierr=NF_REDEF(nid)
776
777	! 2. Define the variable
778	#ifdef NC_DOUBLE
779	ierr = NF_DEF_VAR (nid,adjustl(name),NF_DOUBLE,nbdim,dimids,nvarid)
780	#else
781	ierr = NF_DEF_VAR (nid,adjustl(name),NF_FLOAT,nbdim,dimids,nvarid)
782	#endif
783	if(ierr/=NF_NOERR) then
784	write(,) "def_var_stats: Failed defining variable "//trim(name)
785	write(,) NF_STRERROR(ierr)
786	call abort_physic("def_var_stats","Failed defining "//trim(name),1)
787	endif
788
789	! 3. Write attributes
790	ierr=NF_PUT_ATT_TEXT(nid,nvarid,"title",&
791	len_trim(adjustl(title)),adjustl(title))
792	if(ierr/=NF_NOERR) then
793	write(,) "def_var_stats: Failed writing title attribute for "//trim(name)
794	write(,) NF_STRERROR(ierr)
795	call abort_physic("def_var_stats","Failed writing title for "//trim(name),1)
796	endif
797
798	ierr=NF_PUT_ATT_TEXT(nid,nvarid,"units",&
799	len_trim(adjustl(units)),adjustl(units))
800	if(ierr/=NF_NOERR) then
801	write(,) "def_var_stats: Failed writing units attribute for "//trim(name)
802	write(,) NF_STRERROR(ierr)
803	call abort_physic("def_var_stats","Failed writing units for "//trim(name),1)
804	endif
805
806	! 4. Switch out of NetCDF define mode
807	ierr = NF_ENDDEF(nid)
808
809	end subroutine def_var_stats
810
811	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
812
813	subroutine mkstats(ierr)
814
815	! This routine finalizes tha stats.nc file from sums and sums of squares
816	! to means and standard deviations. The data file is overwritten in place.
817
818	use mod_phys_lmdz_para, only : is_master
819	use mod_grid_phy_lmdz, only : nbp_lon, nbp_lat, nbp_lev, klon_glo
820
821	implicit none
822
823	include "netcdf.inc"
824
825	integer,intent(out) :: ierr ! netCDF return error code
826	integer :: nid,nbvar,i,ndims,lt,nvarid
827	integer, dimension(4) :: id,varid,start,size
828	integer, dimension(5) :: dimids
829	character (len=50) :: name,nameout,units,title
830	real,allocatable :: sum3d(:,:,:),square3d(:,:,:),mean3d(:,:,:),sd3d(:,:,:)
831	real,allocatable :: sum2d(:,:),square2d(:,:),mean2d(:,:),sd2d(:,:)
832	real, dimension(istime) :: time
833	real, dimension(nbp_lat) :: lat
834	real,allocatable :: lon(:)
835	real, dimension(nbp_lev) :: alt
836	logical :: lcopy=.true.
837	!integer :: latid,lonid,altid,timeid
838	integer :: meanid,sdid
839	!integer, dimension(4) :: dimout
840
841	if (is_master) then
842	! only the master needs do all this
843
844	! Incrementation of count for the last step, which is not done in wstats
845	count(istime)=count(istime)+1
846
847	if (klon_glo>1) then
848	allocate(lon(nbp_lon+1))
849	allocate(sum3d(nbp_lon+1,nbp_lat,nbp_lev))
850	allocate(square3d(nbp_lon+1,nbp_lat,nbp_lev))
851	allocate(mean3d(nbp_lon+1,nbp_lat,nbp_lev))
852	allocate(sd3d(nbp_lon+1,nbp_lat,nbp_lev))
853	allocate(sum2d(nbp_lon+1,nbp_lat))
854	allocate(square2d(nbp_lon+1,nbp_lat))
855	allocate(mean2d(nbp_lon+1,nbp_lat))
856	allocate(sd2d(nbp_lon+1,nbp_lat))
857	else ! 1D model case
858	allocate(lon(1))
859	allocate(sum3d(1,1,nbp_lev))
860	allocate(square3d(1,1,nbp_lev))
861	allocate(mean3d(1,1,nbp_lev))
862	allocate(sd3d(1,1,nbp_lev))
863	allocate(sum2d(1,1))
864	allocate(square2d(1,1))
865	allocate(mean2d(1,1))
866	allocate(sd2d(1,1))
867	endif
868
869	ierr = NF_OPEN("stats.nc",NF_WRITE,nid)
870
871	! We catch the id of dimensions of the stats file
872
873	ierr= NF_INQ_DIMID(nid,"latitude",id(1))
874	ierr= NF_INQ_DIMID(nid,"longitude",id(2))
875	ierr= NF_INQ_DIMID(nid,"altitude",id(3))
876	ierr= NF_INQ_DIMID(nid,"Time",id(4))
877
878	ierr= NF_INQ_VARID(nid,"latitude",varid(1))
879	ierr= NF_INQ_VARID(nid,"longitude",varid(2))
880	ierr= NF_INQ_VARID(nid,"altitude",varid(3))
881	ierr= NF_INQ_VARID(nid,"Time",varid(4))
882
883	! Time initialisation
884
885	do i=1,istime
886	time(i)=i*24./istime
887	#ifdef NC_DOUBLE
888	ierr= NF_PUT_VARA_DOUBLE(nid,varid(4),[i],[1],time(i))
889	#else
890	ierr= NF_PUT_VARA_REAL(nid,varid(4),[i],[1],time(i))
891	#endif
892	enddo
893
894	! We catche the values of the variables
895
896	#ifdef NC_DOUBLE
897	ierr = NF_GET_VAR_DOUBLE(nid,varid(1),lat)
898	ierr = NF_GET_VAR_DOUBLE(nid,varid(2),lon)
899	ierr = NF_GET_VAR_DOUBLE(nid,varid(3),alt)
900	#else
901	ierr = NF_GET_VAR_REAL(nid,varid(1),lat)
902	ierr = NF_GET_VAR_REAL(nid,varid(2),lon)
903	ierr = NF_GET_VAR_REAL(nid,varid(3),alt)
904	#endif
905
906	! We catch the number of variables in the stats file
907	ierr = NF_INQ_NVARS(nid,nbvar)
908
909	! to catche the "real" number of variables (without the "additionnal variables")
910	nbvar=(nbvar-4)/2
911
912	do i=1,nbvar
913	nvarid=(i-1)*2+5
914
915	! What's the variable's name?
916	ierr=NF_INQ_VARNAME(nid,nvarid,name)
917	write(,) "OK variable ",name
918	! Its units?
919	units=" "
920	ierr=NF_GET_ATT_TEXT(nid,nvarid,"units",units)
921	! Its title?
922	title=" "
923	ierr=NF_GET_ATT_TEXT(nid,nvarid,"title",title)
924	! Its number of dimensions?
925	ierr=NF_INQ_VARNDIMS(nid,nvarid,ndims)
926	! Its values?
927
928	if(ndims==4) then ! lat, lon, alt & time
929
930	! dimout(1)=lonid
931	! dimout(2)=latid
932	! dimout(3)=altid
933	! dimout(4)=timeid
934
935	if (klon_glo>1) then ! general case
936	size=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
937	else ! 1D model
938	size=(/1,1,nbp_lev,1/)
939	endif
940	do lt=1,istime
941	start=(/1,1,1,lt/)
942	! Extraction of the "source" variables
943	#ifdef NC_DOUBLE
944	ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,size,sum3d)
945	ierr = NF_GET_VARA_DOUBLE(nid,nvarid+1,start,size,square3d)
946	#else
947	ierr = NF_GET_VARA_REAL(nid,nvarid,start,size,sum3d)
948	ierr = NF_GET_VARA_REAL(nid,nvarid+1,start,size,square3d)
949	#endif
950	! Calculation of these nvariables
951	mean3d=sum3d/count(lt)
952	sd3d=sqrt(max(0.,square3d/count(lt)-mean3d**2))
953	! Writing of the variables
954	#ifdef NC_DOUBLE
955	ierr = NF_PUT_VARA_DOUBLE(nid,nvarid,start,size,mean3d)
956	ierr = NF_PUT_VARA_DOUBLE(nid,nvarid+1,start,size,sd3d)
957	#else
958	ierr = NF_PUT_VARA_REAL(nid,nvarid,start,size,mean3d)
959	ierr = NF_PUT_VARA_REAL(nid,nvarid+1,start,size,sd3d)
960	#endif
961	enddo
962
963	else if (ndims.eq.3) then
964
965	! dimout(1)=lonid
966	! dimout(2)=latid
967	! dimout(3)=timeid
968
969	if (klon_glo>1) then ! general case
970	size=(/nbp_lon+1,nbp_lat,1,0/)
971	else
972	size=(/1,1,1,0/)
973	endif
974	do lt=1,istime
975	start=(/1,1,lt,0/)
976	! Extraction of the "source" variables
977	#ifdef NC_DOUBLE
978	ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,size,sum2d)
979	ierr = NF_GET_VARA_DOUBLE(nid,nvarid+1,start,size,square2d)
980	#else
981	ierr = NF_GET_VARA_REAL(nid,nvarid,start,size,sum2d)
982	ierr = NF_GET_VARA_REAL(nid,nvarid+1,start,size,square2d)
983	#endif
984	! Calculation of these variables
985	mean2d=sum2d/count(lt)
986	sd2d=sqrt(max(0.,square2d/count(lt)-mean2d**2))
987	! Writing of the variables
988	#ifdef NC_DOUBLE
989	ierr = NF_PUT_VARA_DOUBLE(nid,nvarid,start,size,mean2d)
990	ierr = NF_PUT_VARA_DOUBLE(nid,nvarid+1,start,size,sd2d)
991	#else
992	ierr = NF_PUT_VARA_REAL(nid,nvarid,start,size,mean2d)
993	ierr = NF_PUT_VARA_REAL(nid,nvarid+1,start,size,sd2d)
994	#endif
995	enddo
996
997	endif
998	enddo
999
1000	ierr= NF_CLOSE(nid)
1001
1002	endif ! of if (is_master)
1003
1004	end subroutine mkstats
1005
1006	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1007
1008	end module wstats_mod

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

Download in other formats:

Original Format