source: LMDZ5/trunk/libf/obsolete/wstats.F90 @ 4002

Last change on this file since 4002 was 2321, checked in by Ehouarn Millour, 9 years ago

Create the "obsolete" directory where old and unused stuff should go. And move some obsolete routines there.
Minor correction in "print_control_mod": use "opened" argument to inquire instead of "exist".
EM
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 8.5 KB
RevLine 
[1486]1subroutine wstats(ngrid,nom,titre,unite,dim,px)
2
3implicit none
4
5#include "dimensions.h"
6#include "statto.h"
7#include "netcdf.inc"
8
9integer,intent(in) :: ngrid
10character (len=*),intent(in) :: nom,titre,unite
11integer,intent(in) :: dim
12real, dimension(ngrid,llm),intent(in) :: px
13integer,parameter :: iip1=iim+1
14integer,parameter :: jjp1=jjm+1
15real, dimension(iip1,jjp1,llm) :: mean3d,sd3d,dx3
16real, dimension(iip1,jjp1) :: mean2d,sd2d,dx2
17character (len=50) :: namebis
18character (len=50), save :: firstvar
19integer :: ierr,varid,nbdim,nid
20integer :: meanid,sdid
21integer, dimension(4)  :: id,start,size
22logical, save :: firstcall=.TRUE.
23integer :: l,i,j,ig0
24integer,save :: index
25
26integer, save :: step=0
27
28
29if (firstcall) then
30   firstcall=.false.
31   firstvar=trim((nom))
32   call inistats(ierr)
33endif
34
35if (firstvar==nom) then ! If we're back to the first variable
36      step = step + 1
37endif
38
39if (mod(step,istats).ne.0) then
40   RETURN
41endif
42
43ierr = NF_OPEN("stats.nc",NF_WRITE,nid)
44
45namebis=trim(nom)
46ierr= NF_INQ_VARID(nid,namebis,meanid)
47
48if (ierr.ne.NF_NOERR) then
49
50   if (firstvar==nom) then
51      index=1
52      count=0
53   endif
54
55
56!declaration de la variable
57
58! choix du nom des coordonnees
59   ierr= NF_INQ_DIMID(nid,"longitude",id(1))
60   ierr= NF_INQ_DIMID(nid,"latitude",id(2))
61   if (dim.eq.3) then
62      ierr= NF_INQ_DIMID(nid,"altitude",id(3))
63      ierr= NF_INQ_DIMID(nid,"Time",id(4))
64      nbdim=4
65   else if (dim==2) then
66      ierr= NF_INQ_DIMID(nid,"Time",id(3))
67      nbdim=3
68   endif
69
70   write (*,*) "====================="
71   write (*,*) "STATS: creation de ",nom
72   namebis=trim(nom)
73   call def_var_stats(nid,namebis,titre,unite,nbdim,id,meanid,ierr)
74   call inivar(nid,meanid,ngrid,dim,index,px,ierr)
75   namebis=trim(nom)//"_sd"
76   call def_var_stats(nid,namebis,trim(titre)//" total standard deviation over the season",unite,nbdim,id,sdid,ierr)
77   call inivar(nid,sdid,ngrid,dim,index,px,ierr)
78
79   ierr= NF_CLOSE(nid)
80   return
81
82else
83   namebis=trim(nom)//"_sd"
84   ierr= NF_INQ_VARID(nid,namebis,sdid)
85
86endif
87
88if (firstvar==nom) then
89   count(index)=count(int(index))+1
90   index=index+1
91   if (index>istime) then
92      index=1
93   endif
94endif
95
96if (count(index)==0) then
97   if (dim.eq.3) then
98      start=(/1,1,1,index/)
99      size=(/iip1,jjp1,llm,1/)
100      mean3d=0
101      sd3d=0
102   else if (dim.eq.2) then
103      start=(/1,1,index,0/)
104      size=(/iip1,jjp1,1,0/)
105      mean2d=0
106      sd2d=0
107   endif
108else
109   if (dim.eq.3) then
110      start=(/1,1,1,index/)
111      size=(/iip1,jjp1,llm,1/)
112#ifdef NC_DOUBLE
113      ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,size,mean3d)
114      ierr = NF_GET_VARA_DOUBLE(nid,sdid,start,size,sd3d)
115#else
116      ierr = NF_GET_VARA_REAL(nid,meanid,start,size,mean3d)
117      ierr = NF_GET_VARA_REAL(nid,sdid,start,size,sd3d)
118#endif
119      if (ierr.ne.NF_NOERR) then
120         write (*,*) NF_STRERROR(ierr)
121         stop ""
122      endif
123
124   else if (dim.eq.2) then
125      start=(/1,1,index,0/)
126      size=(/iip1,jjp1,1,0/)
127#ifdef NC_DOUBLE
128      ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,size,mean2d)
129      ierr = NF_GET_VARA_DOUBLE(nid,sdid,start,size,sd2d)
130#else
131      ierr = NF_GET_VARA_REAL(nid,meanid,start,size,mean2d)
132      ierr = NF_GET_VARA_REAL(nid,sdid,start,size,sd2d)
133#endif
134      if (ierr.ne.NF_NOERR) then
135         write (*,*) NF_STRERROR(ierr)
136         stop ""
137      endif
138   endif
139endif
140
141if (dim.eq.3) then
142
143!  Passage variable physique -->  variable dynamique
144
145   DO l=1,llm
146      DO i=1,iip1
147         dx3(i,1,l)=px(1,l)
148         dx3(i,jjp1,l)=px(ngrid,l)
149      ENDDO
150      DO j=2,jjm
151         ig0= 1+(j-2)*iim
152         DO i=1,iim
153            dx3(i,j,l)=px(ig0+i,l)
154         ENDDO
155         dx3(iip1,j,l)=dx3(1,j,l)
156      ENDDO
157   ENDDO
158
159   mean3d(:,:,:)= mean3d(:,:,:)+dx3(:,:,:)
160   sd3d(:,:,:)= sd3d(:,:,:)+dx3(:,:,:)**2
161
162#ifdef NC_DOUBLE
163   ierr = NF_PUT_VARA_DOUBLE(nid,meanid,start,size,mean3d)
164   ierr = NF_PUT_VARA_DOUBLE(nid,sdid,start,size,sd3d)
165#else
166   ierr = NF_PUT_VARA_REAL(nid,meanid,start,size,mean3d)
167   ierr = NF_PUT_VARA_REAL(nid,sdid,start,size,sd3d)
168#endif
169   if (ierr.ne.NF_NOERR) then
170     write (*,*) NF_STRERROR(ierr)
171     stop ""
172   endif
173
174else if (dim.eq.2) then
175
176!    Passage variable physique -->  physique dynamique
177
178  DO i=1,iip1
179     dx2(i,1)=px(1,1)
180     dx2(i,jjp1)=px(ngrid,1)
181  ENDDO
182  DO j=2,jjm
183     ig0= 1+(j-2)*iim
184     DO i=1,iim
185        dx2(i,j)=px(ig0+i,1)
186     ENDDO
187     dx2(iip1,j)=dx2(1,j)
188  ENDDO
189
190   mean2d(:,:)= mean2d(:,:)+dx2(:,:)
191   sd2d(:,:)= sd2d(:,:)+dx2(:,:)**2
192
193#ifdef NC_DOUBLE
194   ierr = NF_PUT_VARA_DOUBLE(nid,meanid,start,size,mean2d)
195   ierr = NF_PUT_VARA_DOUBLE(nid,sdid,start,size,sd2d)
196#else
197   ierr = NF_PUT_VARA_REAL(nid,meanid,start,size,mean2d)
198   ierr = NF_PUT_VARA_REAL(nid,sdid,start,size,sd2d)
199#endif
200   if (ierr.ne.NF_NOERR) then
201     write (*,*) NF_STRERROR(ierr)
202     stop ""
203   endif
204
205endif
206
207ierr= NF_CLOSE(nid)
208
209end subroutine wstats
210
211!======================================================
212subroutine inivar(nid,varid,ngrid,dim,index,px,ierr)
213
214implicit none
215
216include "dimensions.h"
217!!!!!!include "dimphys.h"
218include "netcdf.inc"
219
220integer, intent(in) :: nid,varid,dim,index,ngrid
221real, dimension(ngrid,llm), intent(in) :: px
222integer, intent(out) :: ierr
223
224integer,parameter :: iip1=iim+1
225integer,parameter :: jjp1=jjm+1
226
227integer :: l,i,j,ig0
228integer, dimension(4) :: start,size
229real, dimension(iip1,jjp1,llm) :: dx3
230real, dimension(iip1,jjp1) :: dx2
231
232if (dim.eq.3) then
233
234   start=(/1,1,1,index/)
235   size=(/iip1,jjp1,llm,1/)
236
237!  Passage variable physique -->  variable dynamique
238
239   DO l=1,llm
240      DO i=1,iip1
241         dx3(i,1,l)=px(1,l)
242         dx3(i,jjp1,l)=px(ngrid,l)
243      ENDDO
244      DO j=2,jjm
245         ig0= 1+(j-2)*iim
246         DO i=1,iim
247            dx3(i,j,l)=px(ig0+i,l)
248         ENDDO
249         dx3(iip1,j,l)=dx3(1,j,l)
250      ENDDO
251   ENDDO
252
253#ifdef NC_DOUBLE
254   ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,size,dx3)
255#else
256   ierr = NF_PUT_VARA_REAL(nid,varid,start,size,dx3)
257#endif
258
259else if (dim.eq.2) then
260
261      start=(/1,1,index,0/)
262      size=(/iip1,jjp1,1,0/)
263
264!    Passage variable physique -->  physique dynamique
265
266  DO i=1,iip1
267     dx2(i,1)=px(1,1)
268     dx2(i,jjp1)=px(ngrid,1)
269  ENDDO
270  DO j=2,jjm
271     ig0= 1+(j-2)*iim
272     DO i=1,iim
273        dx2(i,j)=px(ig0+i,1)
274     ENDDO
275     dx2(iip1,j)=dx2(1,j)
276  ENDDO
277
278#ifdef NC_DOUBLE
279   ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,size,dx2)
280#else
281   ierr = NF_PUT_VARA_REAL(nid,varid,start,size,dx2)
282#endif
283
284endif
285
286end subroutine inivar
287
288!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
289
290subroutine def_var_stats(nid,name,title,units,nbdim,dimids,nvarid,ierr)
291
292! This subroutine defines variable 'name' in a (pre-existing and opened)
293! NetCDF file (known from its NetCDF ID 'nid').
294! The number of dimensions 'nbdim' of the variable, as well as the IDs of
295! corresponding dimensions must be set (in array 'dimids').
296! Upon successfull definition of the variable, 'nvarid' contains the
297! NetCDF ID of the variable.
298! The variables' attributes 'title' (Note that 'long_name' would be more
299! appropriate) and 'units' are also set.
300
301implicit none
302
303#include "netcdf.inc"
304
305integer,intent(in) :: nid ! NetCDF file ID
306character(len=*),intent(in) :: name ! the variable's name
307character(len=*),intent(in) :: title ! 'title' attribute of variable
308character(len=*),intent(in) :: units ! 'units' attribute of variable
309integer,intent(in) :: nbdim ! number of dimensions of the variable
310integer,dimension(nbdim),intent(in) :: dimids ! NetCDF IDs of the dimensions
311                                              ! the variable is defined along
312integer,intent(out) :: nvarid ! NetCDF ID of the variable
313integer,intent(out) :: ierr ! returned NetCDF staus code
314
315! 1. Switch to NetCDF define mode
316ierr=NF_REDEF(nid)
317
318! 2. Define the variable
319#ifdef NC_DOUBLE
320ierr = NF_DEF_VAR (nid,adjustl(name),NF_DOUBLE,nbdim,dimids,nvarid)
321#else
322ierr = NF_DEF_VAR (nid,adjustl(name),NF_FLOAT,nbdim,dimids,nvarid)
323#endif
324if(ierr/=NF_NOERR) then
325   write(*,*) "def_var_stats: Failed defining variable "//trim(name)
326   write(*,*) NF_STRERROR(ierr)
327   stop ""
328endif
329
330! 3. Write attributes
331ierr=NF_PUT_ATT_TEXT(nid,nvarid,"title",&
332                     len_trim(adjustl(title)),adjustl(title))
333if(ierr/=NF_NOERR) then
334   write(*,*) "def_var_stats: Failed writing title attribute for "//trim(name)
335   write(*,*) NF_STRERROR(ierr)
336   stop ""
337endif
338
339ierr=NF_PUT_ATT_TEXT(nid,nvarid,"units",&
340                     len_trim(adjustl(units)),adjustl(units))
341if(ierr/=NF_NOERR) then
342   write(*,*) "def_var_stats: Failed writing units attribute for "//trim(name)
343   write(*,*) NF_STRERROR(ierr)
344   stop ""
345endif
346
347! 4. Switch out of NetCDF define mode
348ierr = NF_ENDDEF(nid)
349
350end subroutine def_var_stats
351
Note: See TracBrowser for help on using the repository browser.