source: trunk/LMDZ.MARS/libf/phymars/wstats.F90 @ 804

Last change on this file since 804 was 719, checked in by tnavarro, 13 years ago

use of stats.def, same as diagfi.def, and some cleaning of outputs in physiq.F
File size: 9.9 KB
RevLine 
[38]1subroutine wstats(ngrid,nom,titre,unite,dim,px)
2
3implicit none
4
5#include "dimensions.h"
6#include "dimphys.h"
7#include "statto.h"
8#include "netcdf.inc"
9
10integer,intent(in) :: ngrid
11character (len=*),intent(in) :: nom,titre,unite
12integer,intent(in) :: dim
13real, dimension(ngrid,llm),intent(in) :: px
14integer,parameter :: iip1=iim+1
15integer,parameter :: jjp1=jjm+1
16real, dimension(iip1,jjp1,llm) :: mean3d,sd3d,dx3
17real, dimension(iip1,jjp1) :: mean2d,sd2d,dx2
18character (len=50) :: namebis
19character (len=50), save :: firstvar
[719]20integer :: ierr,ierr2,varid,nbdim,nid
[38]21integer :: meanid,sdid
22integer, dimension(4)  :: id,start,size
23logical, save :: firstcall=.TRUE.
24integer :: l,i,j,ig0
25integer,save :: index
[719]26! Added to use stats.def to select output variable
27logical,save :: stats_def
28logical :: getout
29integer,save :: n_nom_def
30integer :: n
31integer,parameter :: n_nom_def_max=99
32character(len=20),save :: nom_def(n_nom_def_max)
33logical, save :: notfoundfirst=.TRUE.
[38]34
35integer, save :: step=0
[719]36     
37     
[38]38
39if (firstcall) then
40   firstcall=.false.
41   call inistats(ierr)
[719]42! at very first call, check if there is a "stats.def" to use and read it
43   open(99,file="stats.def",status='old',form='formatted',iostat=ierr2)
44   if(ierr2.eq.0) then
45     stats_def=.true.
46     write(*,*) "******************"
47     write(*,*) "Reading stats.def"
48     write(*,*) "******************"
49     do n=1,n_nom_def_max
50          read(99,fmt='(a)',end=88) nom_def(n)
51          write(*,*) 'Output in stats: ', trim(nom_def(n))
52     end do
5388   continue
54     if (n.ge.n_nom_def_max) then
55        write(*,*)"n_nom_def_max too small in wstats.F90:",n
56        stop
57     end if
58     n_nom_def=n-1
59     close(99)
60   else
61     firstvar=trim((nom))
62     stats_def=.false.
63     notfoundfirst=.false.
64   endif
[38]65endif
66
[719]67! find firstvar that is in stats.def
68if (notfoundfirst) then
69   do n=1,n_nom_def_max
70     if(trim((nom_def(n))).eq.trim((nom))) then
71        firstvar=trim((nom))
72        notfoundfirst=.false.
73     endif
74   end do
75endif
76
77! get out of wstats if there is stats.def AND variable not listed
78if (stats_def) then
79   getout=.true.
80   do n=1,n_nom_def
81       if(trim(nom_def(n)).eq.nom) getout=.false.
82   end do
83   if (getout) return
84end if
85
[38]86if (firstvar==nom) then ! If we're back to the first variable
87      step = step + 1
88endif
89
90if (mod(step,istats).ne.0) then
91   RETURN
92endif
93
94ierr = NF_OPEN("stats.nc",NF_WRITE,nid)
95
96namebis=trim(nom)
97ierr= NF_INQ_VARID(nid,namebis,meanid)
98
99if (ierr.ne.NF_NOERR) then
100
101   if (firstvar==nom) then
102      index=1
103      count=0
104   endif
105
106
107!declaration de la variable
108
109! choix du nom des coordonnees
110   ierr= NF_INQ_DIMID(nid,"longitude",id(1))
111   ierr= NF_INQ_DIMID(nid,"latitude",id(2))
112   if (dim.eq.3) then
113      ierr= NF_INQ_DIMID(nid,"altitude",id(3))
114      ierr= NF_INQ_DIMID(nid,"Time",id(4))
115      nbdim=4
116   else if (dim==2) then
117      ierr= NF_INQ_DIMID(nid,"Time",id(3))
118      nbdim=3
119   endif
120
121   write (*,*) "====================="
122   write (*,*) "STATS: creation de ",nom
123   namebis=trim(nom)
124   call def_var_stats(nid,namebis,titre,unite,nbdim,id,meanid,ierr)
125   call inivar(nid,meanid,ngrid,dim,index,px,ierr)
126   namebis=trim(nom)//"_sd"
127   call def_var_stats(nid,namebis,trim(titre)//" total standard deviation over the season",unite,nbdim,id,sdid,ierr)
128   call inivar(nid,sdid,ngrid,dim,index,px,ierr)
129
130   ierr= NF_CLOSE(nid)
131   return
132
133else
134   namebis=trim(nom)//"_sd"
135   ierr= NF_INQ_VARID(nid,namebis,sdid)
136
137endif
138
139if (firstvar==nom) then
140   count(index)=count(int(index))+1
141   index=index+1
142   if (index>istime) then
143      index=1
144   endif
145endif
146
147if (count(index)==0) then
148   if (dim.eq.3) then
149      start=(/1,1,1,index/)
150      size=(/iip1,jjp1,llm,1/)
151      mean3d=0
152      sd3d=0
153   else if (dim.eq.2) then
154      start=(/1,1,index,0/)
155      size=(/iip1,jjp1,1,0/)
156      mean2d=0
157      sd2d=0
158   endif
159else
160   if (dim.eq.3) then
161      start=(/1,1,1,index/)
162      size=(/iip1,jjp1,llm,1/)
163#ifdef NC_DOUBLE
164      ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,size,mean3d)
165      ierr = NF_GET_VARA_DOUBLE(nid,sdid,start,size,sd3d)
166#else
167      ierr = NF_GET_VARA_REAL(nid,meanid,start,size,mean3d)
168      ierr = NF_GET_VARA_REAL(nid,sdid,start,size,sd3d)
169#endif
170      if (ierr.ne.NF_NOERR) then
171         write (*,*) NF_STRERROR(ierr)
172         stop ""
173      endif
174
175   else if (dim.eq.2) then
176      start=(/1,1,index,0/)
177      size=(/iip1,jjp1,1,0/)
178#ifdef NC_DOUBLE
179      ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,size,mean2d)
180      ierr = NF_GET_VARA_DOUBLE(nid,sdid,start,size,sd2d)
181#else
182      ierr = NF_GET_VARA_REAL(nid,meanid,start,size,mean2d)
183      ierr = NF_GET_VARA_REAL(nid,sdid,start,size,sd2d)
184#endif
185      if (ierr.ne.NF_NOERR) then
186         write (*,*) NF_STRERROR(ierr)
187         stop ""
188      endif
189   endif
190endif
191
192if (dim.eq.3) then
193
194!  Passage variable physique -->  variable dynamique
195
196   DO l=1,llm
197      DO i=1,iip1
198         dx3(i,1,l)=px(1,l)
199         dx3(i,jjp1,l)=px(ngrid,l)
200      ENDDO
201      DO j=2,jjm
202         ig0= 1+(j-2)*iim
203         DO i=1,iim
204            dx3(i,j,l)=px(ig0+i,l)
205         ENDDO
206         dx3(iip1,j,l)=dx3(1,j,l)
207      ENDDO
208   ENDDO
209
210   mean3d(:,:,:)= mean3d(:,:,:)+dx3(:,:,:)
211   sd3d(:,:,:)= sd3d(:,:,:)+dx3(:,:,:)**2
212
213#ifdef NC_DOUBLE
214   ierr = NF_PUT_VARA_DOUBLE(nid,meanid,start,size,mean3d)
215   ierr = NF_PUT_VARA_DOUBLE(nid,sdid,start,size,sd3d)
216#else
217   ierr = NF_PUT_VARA_REAL(nid,meanid,start,size,mean3d)
218   ierr = NF_PUT_VARA_REAL(nid,sdid,start,size,sd3d)
219#endif
220   if (ierr.ne.NF_NOERR) then
221     write (*,*) NF_STRERROR(ierr)
222     stop ""
223   endif
224
225else if (dim.eq.2) then
226
227!    Passage variable physique -->  physique dynamique
228
229  DO i=1,iip1
230     dx2(i,1)=px(1,1)
231     dx2(i,jjp1)=px(ngrid,1)
232  ENDDO
233  DO j=2,jjm
234     ig0= 1+(j-2)*iim
235     DO i=1,iim
236        dx2(i,j)=px(ig0+i,1)
237     ENDDO
238     dx2(iip1,j)=dx2(1,j)
239  ENDDO
240
241   mean2d(:,:)= mean2d(:,:)+dx2(:,:)
242   sd2d(:,:)= sd2d(:,:)+dx2(:,:)**2
243
244#ifdef NC_DOUBLE
245   ierr = NF_PUT_VARA_DOUBLE(nid,meanid,start,size,mean2d)
246   ierr = NF_PUT_VARA_DOUBLE(nid,sdid,start,size,sd2d)
247#else
248   ierr = NF_PUT_VARA_REAL(nid,meanid,start,size,mean2d)
249   ierr = NF_PUT_VARA_REAL(nid,sdid,start,size,sd2d)
250#endif
251   if (ierr.ne.NF_NOERR) then
252     write (*,*) NF_STRERROR(ierr)
253     stop ""
254   endif
255
256endif
257
258ierr= NF_CLOSE(nid)
259
260end subroutine wstats
261
262!======================================================
263subroutine inivar(nid,varid,ngrid,dim,index,px,ierr)
264
265implicit none
266
267include "dimensions.h"
268include "dimphys.h"
269include "netcdf.inc"
270
271integer, intent(in) :: nid,varid,dim,index,ngrid
272real, dimension(ngrid,llm), intent(in) :: px
273integer, intent(out) :: ierr
274
275integer,parameter :: iip1=iim+1
276integer,parameter :: jjp1=jjm+1
277
278integer :: l,i,j,ig0
279integer, dimension(4) :: start,size
280real, dimension(iip1,jjp1,llm) :: dx3
281real, dimension(iip1,jjp1) :: dx2
282
283if (dim.eq.3) then
284
285   start=(/1,1,1,index/)
286   size=(/iip1,jjp1,llm,1/)
287
288!  Passage variable physique -->  variable dynamique
289
290   DO l=1,llm
291      DO i=1,iip1
292         dx3(i,1,l)=px(1,l)
293         dx3(i,jjp1,l)=px(ngrid,l)
294      ENDDO
295      DO j=2,jjm
296         ig0= 1+(j-2)*iim
297         DO i=1,iim
298            dx3(i,j,l)=px(ig0+i,l)
299         ENDDO
300         dx3(iip1,j,l)=dx3(1,j,l)
301      ENDDO
302   ENDDO
303
304#ifdef NC_DOUBLE
305   ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,size,dx3)
306#else
307   ierr = NF_PUT_VARA_REAL(nid,varid,start,size,dx3)
308#endif
309
310else if (dim.eq.2) then
311
312      start=(/1,1,index,0/)
313      size=(/iip1,jjp1,1,0/)
314
315!    Passage variable physique -->  physique dynamique
316
317  DO i=1,iip1
318     dx2(i,1)=px(1,1)
319     dx2(i,jjp1)=px(ngrid,1)
320  ENDDO
321  DO j=2,jjm
322     ig0= 1+(j-2)*iim
323     DO i=1,iim
324        dx2(i,j)=px(ig0+i,1)
325     ENDDO
326     dx2(iip1,j)=dx2(1,j)
327  ENDDO
328
329#ifdef NC_DOUBLE
330   ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,size,dx2)
331#else
332   ierr = NF_PUT_VARA_REAL(nid,varid,start,size,dx2)
333#endif
334
335endif
336
337end subroutine inivar
338
339!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
340
341subroutine def_var_stats(nid,name,title,units,nbdim,dimids,nvarid,ierr)
342
343! This subroutine defines variable 'name' in a (pre-existing and opened)
344! NetCDF file (known from its NetCDF ID 'nid').
345! The number of dimensions 'nbdim' of the variable, as well as the IDs of
346! corresponding dimensions must be set (in array 'dimids').
347! Upon successfull definition of the variable, 'nvarid' contains the
348! NetCDF ID of the variable.
349! The variables' attributes 'title' (Note that 'long_name' would be more
350! appropriate) and 'units' are also set.
351
352implicit none
353
354#include "netcdf.inc"
355
356integer,intent(in) :: nid ! NetCDF file ID
357character(len=*),intent(in) :: name ! the variable's name
358character(len=*),intent(in) :: title ! 'title' attribute of variable
359character(len=*),intent(in) :: units ! 'units' attribute of variable
360integer,intent(in) :: nbdim ! number of dimensions of the variable
361integer,dimension(nbdim),intent(in) :: dimids ! NetCDF IDs of the dimensions
362                                              ! the variable is defined along
363integer,intent(out) :: nvarid ! NetCDF ID of the variable
364integer,intent(out) :: ierr ! returned NetCDF staus code
365
366! 1. Switch to NetCDF define mode
367ierr=NF_REDEF(nid)
368
369! 2. Define the variable
370#ifdef NC_DOUBLE
371ierr = NF_DEF_VAR (nid,adjustl(name),NF_DOUBLE,nbdim,dimids,nvarid)
372#else
373ierr = NF_DEF_VAR (nid,adjustl(name),NF_FLOAT,nbdim,dimids,nvarid)
374#endif
375if(ierr/=NF_NOERR) then
376   write(*,*) "def_var_stats: Failed defining variable "//trim(name)
377   write(*,*) NF_STRERROR(ierr)
378   stop ""
379endif
380
381! 3. Write attributes
382ierr=NF_PUT_ATT_TEXT(nid,nvarid,"title",&
383                     len_trim(adjustl(title)),adjustl(title))
384if(ierr/=NF_NOERR) then
385   write(*,*) "def_var_stats: Failed writing title attribute for "//trim(name)
386   write(*,*) NF_STRERROR(ierr)
387   stop ""
388endif
389
390ierr=NF_PUT_ATT_TEXT(nid,nvarid,"units",&
391                     len_trim(adjustl(units)),adjustl(units))
392if(ierr/=NF_NOERR) then
393   write(*,*) "def_var_stats: Failed writing units attribute for "//trim(name)
394   write(*,*) NF_STRERROR(ierr)
395   stop ""
396endif
397
398! 4. Switch out of NetCDF define mode
399ierr = NF_ENDDEF(nid)
400
401end subroutine def_var_stats
402
Note: See TracBrowser for help on using the repository browser.