Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

wstats_mod.F90

Last change on this file 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

Rev	Line
[2559]	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
[2573]	20	!$OMP THREADPRIVATE(count)
[2559]	21
	22	contains
	23
[38]	24	subroutine wstats(ngrid,nom,titre,unite,dim,px)
	25
[2559]	26	! Main routine to call to trigger writing a given field to the "stats.nc" file
	27
[3282]	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
[1528]	30	use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo, &
[2563]	31	nbp_lon, nbp_lat, nbp_lev, &
	32	grid_type, unstructured
[38]	33	implicit none
	34
[1528]	35	include "netcdf.inc"
[38]	36
	37	integer,intent(in) :: ngrid
	38	character (len=*),intent(in) :: nom,titre,unite
	39	integer,intent(in) :: dim
[1528]	40	real,intent(in) :: px(ngrid,nbp_lev)
[1532]	41	real,allocatable,save :: mean3d(:,:,:),sd3d(:,:,:),dx3(:,:,:)
	42	real,allocatable,save :: mean2d(:,:),sd2d(:,:),dx2(:,:)
[38]	43	character (len=50) :: namebis
	44	character (len=50), save :: firstvar
[2616]	45	!$OMP THREADPRIVATE(mean3d,sd3d,dx3,mean2d,sd2d,dx2,firstvar)
[1130]	46	integer :: ierr,varid,nbdim,nid
[38]	47	integer :: meanid,sdid
[1130]	48	integer, dimension(4) :: id,start,sizes
[38]	49	logical, save :: firstcall=.TRUE.
[1130]	50	integer,save :: indx
[38]	51	integer, save :: step=0
[1532]	52	!$OMP THREADPRIVATE(firstcall,indx,step)
[2563]	53	integer :: l,i,j,ig0,n
[38]	54
[2563]	55	! Added to read an optional stats.def file to select outputs
[3282]	56	logical,save :: check_stats_def=.TRUE. ! to keep track of if stats.def was checked
[2563]	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
[3282]	63	!$OMP THREADPRIVATE(check_stats_def,stats_def,n_name_def,name_def)
[2616]	64
[1130]	65	! Added to work in parallel mode
	66	#ifdef CPP_PARA
[1528]	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
[1130]	69	real px2_glop(klon_glo) ! to store a 2D data set on global physics grid
[1528]	70	real px2_glo(nbp_lon,nbp_lat) ! to store a 2D data set on global lonxlat grid
[1130]	71	real px2(ngrid)
[1528]	72	real px3(ngrid,nbp_lev)
[1130]	73	#else
	74	! When not running in parallel mode:
[1528]	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
[1130]	77	real px2_glop(ngrid) ! to store a 2D data set on global physics grid
[1528]	78	real px2_glo(nbp_lon,nbp_lat) ! to store a 2D data set on global lonxlat grid
[1130]	79	#endif
	80
[2563]	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
[3282]	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
[2563]	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)
[3282]	118	else
[2563]	119	stats_def=.false. ! no stats.def file; output all fields sent to wstats
[3282]	120	endif ! of if (ierr.eq.0)
[2563]	121
[3282]	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
[38]	166	call inistats(ierr)
[3282]	167
	168	! allocate arrays that will be needed later
[1532]	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
[3282]	184
[2563]	185	endif ! of if (firstcall)
[38]	186
[1130]	187	if (firstvar==nom) then ! If we're back to the first variable, increment time counter
[38]	188	step = step + 1
	189	endif
	190
	191	if (mod(step,istats).ne.0) then
[1130]	192	! if its not time to write to file, exit
[38]	193	RETURN
	194	endif
	195
[2563]	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
[1130]	213	! collect fields on a global physics grid
	214	#ifdef CPP_PARA
	215	if (dim.eq.3) then
[1528]	216	px3(1:ngrid,1:nbp_lev)=px(1:ngrid,1:nbp_lev)
[1130]	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)
[2997]	222	if (klon_glo>1) then ! general case (unless in 1D)
[1130]	223	! copy dx3_glo() to dx3(:) and add redundant longitude
[1528]	224	dx3(1:nbp_lon,:,:)=px3_glo(1:nbp_lon,:,:)
	225	dx3(nbp_lon+1,:,:)=dx3(1,:,:)
[2997]	226	endif
[1130]	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
[2997]	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
[1130]	243	!$OMP END MASTER
	244	!$OMP BARRIER
	245	endif
	246	#else
	247	if (dim.eq.3) then
[1528]	248	px3_glop(:,1:nbp_lev)=px(:,1:nbp_lev)
[1130]	249	! Passage variable physique --> variable dynamique
[1528]	250	DO l=1,nbp_lev
	251	DO i=1,nbp_lon
[1130]	252	px3_glo(i,1,l)=px(1,l)
[1528]	253	px3_glo(i,nbp_lat,l)=px(ngrid,l)
[1130]	254	ENDDO
[1528]	255	DO j=2,nbp_lat-1
	256	ig0= 1+(j-2)*nbp_lon
	257	DO i=1,nbp_lon
[1130]	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
[1528]	265	DO i=1,nbp_lon
[1130]	266	px2_glo(i,1)=px(1,1)
[1528]	267	px2_glo(i,nbp_lat)=px(ngrid,1)
[1130]	268	ENDDO
[1528]	269	DO j=2,nbp_lat-1
	270	ig0= 1+(j-2)*nbp_lon
	271	DO i=1,nbp_lon
[1130]	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
[38]	283	ierr = NF_OPEN("stats.nc",NF_WRITE,nid)
	284
	285	namebis=trim(nom)
[1130]	286
	287	! test: check if that variable already exists in file
[38]	288	ierr= NF_INQ_VARID(nid,namebis,meanid)
	289
	290	if (ierr.ne.NF_NOERR) then
[1130]	291	! variable not in file, create/define it
[38]	292	if (firstvar==nom) then
[1130]	293	indx=1
	294	count(:)=0
[38]	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 (,) "====================="
[2311]	313	write (,) "STATS: creating ",nom
[38]	314	namebis=trim(nom)
[1560]	315	call def_var_stats(nid,namebis,titre,unite,nbdim,id,meanid,ierr)
[1130]	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
[38]	321	namebis=trim(nom)//"_sd"
[1560]	322	call def_var_stats(nid,namebis,trim(titre)//" total standard deviation over the season",unite,nbdim,id,sdid,ierr)
[1130]	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
[38]	328
	329	ierr= NF_CLOSE(nid)
	330	return
	331
	332	else
[1130]	333	! variable found in file
[38]	334	namebis=trim(nom)//"_sd"
	335	ierr= NF_INQ_VARID(nid,namebis,sdid)
	336
	337	endif
	338
	339	if (firstvar==nom) then
[1130]	340	count(indx)=count(int(indx))+1
	341	indx=indx+1
	342	if (indx>istime) then
	343	indx=1
[38]	344	endif
	345	endif
	346
[1130]	347	if (count(indx)==0) then
	348	! very first time we write the variable to file
[38]	349	if (dim.eq.3) then
[1130]	350	start=(/1,1,1,indx/)
[1532]	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
[1130]	356	mean3d(:,:,:)=0
	357	sd3d(:,:,:)=0
[38]	358	else if (dim.eq.2) then
[1130]	359	start=(/1,1,indx,0/)
[1532]	360	if (klon_glo>1) then !general case
[1689]	361	sizes=(/nbp_lon+1,nbp_lat,1,0/)
[1532]	362	else
	363	sizes=(/1,1,1,0/)
	364	endif
[1130]	365	mean2d(:,:)=0
	366	sd2d(:,:)=0
[38]	367	endif
	368	else
[1130]	369	! load values from file
[38]	370	if (dim.eq.3) then
[1130]	371	start=(/1,1,1,indx/)
[1532]	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
[38]	377	#ifdef NC_DOUBLE
[1130]	378	ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,sizes,mean3d)
	379	ierr = NF_GET_VARA_DOUBLE(nid,sdid,start,sizes,sd3d)
[38]	380	#else
[1130]	381	ierr = NF_GET_VARA_REAL(nid,meanid,start,sizes,mean3d)
	382	ierr = NF_GET_VARA_REAL(nid,sdid,start,sizes,sd3d)
[38]	383	#endif
	384	if (ierr.ne.NF_NOERR) then
[1532]	385	write (,) "wstats error reading :",trim(nom)
[38]	386	write (,) NF_STRERROR(ierr)
[2311]	387	call abort_physic("wstats","Failed reading "//trim(nom),1)
[38]	388	endif
	389
	390	else if (dim.eq.2) then
[1130]	391	start=(/1,1,indx,0/)
[1532]	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
[38]	397	#ifdef NC_DOUBLE
[1130]	398	ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,sizes,mean2d)
	399	ierr = NF_GET_VARA_DOUBLE(nid,sdid,start,sizes,sd2d)
[38]	400	#else
[1130]	401	ierr = NF_GET_VARA_REAL(nid,meanid,start,sizes,mean2d)
	402	ierr = NF_GET_VARA_REAL(nid,sdid,start,sizes,sd2d)
[38]	403	#endif
	404	if (ierr.ne.NF_NOERR) then
[1532]	405	write (,) "wstats error reading :",trim(nom)
[38]	406	write (,) NF_STRERROR(ierr)
[2311]	407	call abort_physic("wstats","Failed reading "//trim(nom),1)
[38]	408	endif
	409	endif
[1130]	410	endif ! of if (count(indx)==0)
[38]	411
[1130]	412
	413	! 2.1. Build dx* (data on lon-lat grid, with redundant longitude)
	414
[38]	415	if (dim.eq.3) then
[1528]	416	dx3(1:nbp_lon,:,:)=px3_glo(:,:,:)
[1532]	417	IF (klon_glo>1) THEN ! in 3D, add redundant longitude point
	418	dx3(nbp_lon+1,:,:)=dx3(1,:,:)
	419	ENDIF
[1130]	420	else ! dim.eq.2
[1528]	421	dx2(1:nbp_lon,:)=px2_glo(:,:)
[1532]	422	IF (klon_glo>1) THEN ! in 3D, add redundant longitude point
	423	dx2(nbp_lon+1,:)=dx2(1,:)
	424	ENDIF
[1130]	425	endif
[38]	426
	427
[1130]	428	! 2.2. Add current values to previously stored sums
[38]	429
[1130]	430	if (dim.eq.3) then
[38]	431
[1130]	432	mean3d(:,:,:)=mean3d(:,:,:)+dx3(:,:,:)
	433	sd3d(:,:,:)=sd3d(:,:,:)+dx3(:,:,:)**2
	434
[38]	435	#ifdef NC_DOUBLE
[1130]	436	ierr = NF_PUT_VARA_DOUBLE(nid,meanid,start,sizes,mean3d)
	437	ierr = NF_PUT_VARA_DOUBLE(nid,sdid,start,sizes,sd3d)
[38]	438	#else
[1130]	439	ierr = NF_PUT_VARA_REAL(nid,meanid,start,sizes,mean3d)
	440	ierr = NF_PUT_VARA_REAL(nid,sdid,start,sizes,sd3d)
[38]	441	#endif
[1532]	442	if (ierr.ne.NF_NOERR) then
	443	write (,) "wstats error writing :",trim(nom)
	444	write (,) NF_STRERROR(ierr)
[2311]	445	call abort_physic("wstats","Failed writing "//trim(nom),1)
[1532]	446	endif
[38]	447
	448	else if (dim.eq.2) then
	449
	450	mean2d(:,:)= mean2d(:,:)+dx2(:,:)
[1130]	451	sd2d(:,:)=sd2d(:,:)+dx2(:,:)**2
[38]	452
	453	#ifdef NC_DOUBLE
[1130]	454	ierr = NF_PUT_VARA_DOUBLE(nid,meanid,start,sizes,mean2d)
	455	ierr = NF_PUT_VARA_DOUBLE(nid,sdid,start,sizes,sd2d)
[38]	456	#else
[1130]	457	ierr = NF_PUT_VARA_REAL(nid,meanid,start,sizes,mean2d)
	458	ierr = NF_PUT_VARA_REAL(nid,sdid,start,sizes,sd2d)
[38]	459	#endif
[1532]	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)
[2311]	467	call abort_physic("wstats","Failed writing "//trim(nom),1)
[1532]	468	endif
[38]	469
[1130]	470	endif ! of if (dim.eq.3) elseif (dim.eq.2)
[38]	471
[1130]	472	ierr= NF_CLOSE(nid)
	473	endif ! of if (is_master)
[38]	474
	475	end subroutine wstats
	476
[2559]	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
[2573]	591	ierr = NF_DEF_VAR (nid,"altitude", NF_DOUBLE, 1,[idim_llm],nvarid)
[2559]	592	#else
[2573]	593	ierr = NF_DEF_VAR (nid,"altitude", NF_FLOAT, 1,[idim_llm],nvarid)
[2559]	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
[1130]	628	subroutine inivar(nid,varid,ngrid,dim,indx,px,ierr)
[2559]	629
	630	! routine to initialize writing a variable in the stats file
	631
[1532]	632	use mod_grid_phy_lmdz, only : nbp_lon, nbp_lat, nbp_lev, klon_glo
[38]	633
	634	implicit none
	635
	636	include "netcdf.inc"
	637
[1130]	638	integer, intent(in) :: nid,varid,dim,indx,ngrid
[1528]	639	real, dimension(ngrid,nbp_lev), intent(in) :: px
[38]	640	integer, intent(out) :: ierr
	641
	642	integer :: l,i,j,ig0
[1130]	643	integer, dimension(4) :: start,sizes
[1528]	644	real, dimension(nbp_lon+1,nbp_lat,nbp_lev) :: dx3
	645	real, dimension(nbp_lon+1,nbp_lat) :: dx2
[1532]	646	real :: dx3_1d(nbp_lev) ! for 1D outputs
	647	real :: dx2_1d ! for 1D outputs
[38]	648
	649	if (dim.eq.3) then
	650
[1130]	651	start=(/1,1,1,indx/)
[1532]	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
[38]	657
	658	! Passage variable physique --> variable dynamique
	659
[1532]	660	if (klon_glo>1) then ! general case
	661	DO l=1,nbp_lev
[1528]	662	DO i=1,nbp_lon+1
[38]	663	dx3(i,1,l)=px(1,l)
[1528]	664	dx3(i,nbp_lat,l)=px(ngrid,l)
[38]	665	ENDDO
[1528]	666	DO j=2,nbp_lat-1
	667	ig0= 1+(j-2)*nbp_lon
	668	DO i=1,nbp_lon
[38]	669	dx3(i,j,l)=px(ig0+i,l)
	670	ENDDO
[1528]	671	dx3(nbp_lon+1,j,l)=dx3(1,j,l)
[38]	672	ENDDO
[1532]	673	ENDDO
	674	else ! 1D model case
	675	dx3_1d(1:nbp_lev)=px(1,1:nbp_lev)
	676	endif
[38]	677
	678	#ifdef NC_DOUBLE
[1532]	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
[38]	684	#else
[1532]	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
[38]	690	#endif
[1532]	691	if (ierr.ne.NF_NOERR) then
	692	write (,) "inivar error writing variable"
	693	write (,) NF_STRERROR(ierr)
[2311]	694	call abort_physic("inivar","error writing variable",1)
[1532]	695	endif
[38]	696
	697	else if (dim.eq.2) then
	698
[1130]	699	start=(/1,1,indx,0/)
[1532]	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
[38]	705
	706	! Passage variable physique --> physique dynamique
	707
[1532]	708	if (klon_glo>1) then ! general case
[1528]	709	DO i=1,nbp_lon+1
[38]	710	dx2(i,1)=px(1,1)
[1528]	711	dx2(i,nbp_lat)=px(ngrid,1)
[38]	712	ENDDO
[1528]	713	DO j=2,nbp_lat-1
	714	ig0= 1+(j-2)*nbp_lon
	715	DO i=1,nbp_lon
[38]	716	dx2(i,j)=px(ig0+i,1)
	717	ENDDO
[1528]	718	dx2(nbp_lon+1,j)=dx2(1,j)
[38]	719	ENDDO
[1532]	720	else ! 1D model case
	721	dx2_1d=px(1,1)
	722	endif
	723
[38]	724	#ifdef NC_DOUBLE
[1532]	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
[38]	730	#else
[1532]	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
[38]	736	#endif
[1532]	737	if (ierr.ne.NF_NOERR) then
	738	write (,) "inivar error writing variable"
	739	write (,) NF_STRERROR(ierr)
[2311]	740	call abort_physic("inivar","error writing variable",1)
[1532]	741	endif
[38]	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
[1528]	762	include "netcdf.inc"
[38]	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)
[2311]	786	call abort_physic("def_var_stats","Failed defining "//trim(name),1)
[38]	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)
[2311]	795	call abort_physic("def_var_stats","Failed writing title for "//trim(name),1)
[38]	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)
[2311]	803	call abort_physic("def_var_stats","Failed writing units for "//trim(name),1)
[38]	804	endif
	805
	806	! 4. Switch out of NetCDF define mode
	807	ierr = NF_ENDDEF(nid)
	808
	809	end subroutine def_var_stats
	810
[2559]	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
[2573]	888	ierr= NF_PUT_VARA_DOUBLE(nid,varid(4),[i],[1],time(i))
[2559]	889	#else
[2573]	890	ierr= NF_PUT_VARA_REAL(nid,varid(4),[i],[1],time(i))
[2559]	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
[2573]	913	nvarid=(i-1)*2+5
[2559]	914
	915	! What's the variable's name?
[2573]	916	ierr=NF_INQ_VARNAME(nid,nvarid,name)
[2559]	917	write(,) "OK variable ",name
	918	! Its units?
	919	units=" "
[2573]	920	ierr=NF_GET_ATT_TEXT(nid,nvarid,"units",units)
[2559]	921	! Its title?
	922	title=" "
[2573]	923	ierr=NF_GET_ATT_TEXT(nid,nvarid,"title",title)
[2559]	924	! Its number of dimensions?
[2573]	925	ierr=NF_INQ_VARNDIMS(nid,nvarid,ndims)
[2559]	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
[2573]	944	ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,size,sum3d)
	945	ierr = NF_GET_VARA_DOUBLE(nid,nvarid+1,start,size,square3d)
[2559]	946	#else
[2573]	947	ierr = NF_GET_VARA_REAL(nid,nvarid,start,size,sum3d)
	948	ierr = NF_GET_VARA_REAL(nid,nvarid+1,start,size,square3d)
[2559]	949	#endif
[2573]	950	! Calculation of these nvariables
[2559]	951	mean3d=sum3d/count(lt)
	952	sd3d=sqrt(max(0.,square3d/count(lt)-mean3d**2))
	953	! Writing of the variables
	954	#ifdef NC_DOUBLE
[2573]	955	ierr = NF_PUT_VARA_DOUBLE(nid,nvarid,start,size,mean3d)
	956	ierr = NF_PUT_VARA_DOUBLE(nid,nvarid+1,start,size,sd3d)
[2559]	957	#else
[2573]	958	ierr = NF_PUT_VARA_REAL(nid,nvarid,start,size,mean3d)
	959	ierr = NF_PUT_VARA_REAL(nid,nvarid+1,start,size,sd3d)
[2559]	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
[2573]	978	ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,size,sum2d)
	979	ierr = NF_GET_VARA_DOUBLE(nid,nvarid+1,start,size,square2d)
[2559]	980	#else
[2573]	981	ierr = NF_GET_VARA_REAL(nid,nvarid,start,size,sum2d)
	982	ierr = NF_GET_VARA_REAL(nid,nvarid+1,start,size,square2d)
[2559]	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
[2573]	989	ierr = NF_PUT_VARA_DOUBLE(nid,nvarid,start,size,mean2d)
	990	ierr = NF_PUT_VARA_DOUBLE(nid,nvarid+1,start,size,sd2d)
[2559]	991	#else
[2573]	992	ierr = NF_PUT_VARA_REAL(nid,nvarid,start,size,mean2d)
	993	ierr = NF_PUT_VARA_REAL(nid,nvarid+1,start,size,sd2d)
[2559]	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