source: trunk/LMDZ.GENERIC/libf/phystd/wstats.F90 @ 1309

Last change on this file since 1309 was 1308, checked in by emillour, 10 years ago

Generic GCM:
Some cleanup to simplify dynamics/physics interactions by getting rid
of dimphys.h (i.e. the nlayermx parameter) and minimizing use of
dimension.h in the physics.
EM
File size: 11.1 KB
Line 
1subroutine wstats(ngrid,nom,titre,unite,dim,px)
2
3use mod_phys_lmdz_para, only : is_mpi_root, is_master, gather, klon_mpi_begin
4use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
5
6implicit none
7
8#include "dimensions.h"
9!#include "dimphys.h"
10#include "comconst.h"
11#include "statto.h"
12#include "netcdf.inc"
13
14integer,intent(in) :: ngrid
15character (len=*),intent(in) :: nom,titre,unite
16integer,intent(in) :: dim
17integer,parameter :: iip1=iim+1
18integer,parameter :: jjp1=jjm+1
19real,intent(in) :: px(ngrid,llm)
20real, dimension(iip1,jjp1,llm) :: mean3d,sd3d,dx3
21real, dimension(iip1,jjp1) :: mean2d,sd2d,dx2
22character (len=50) :: namebis
23character (len=50), save :: firstvar
24integer :: ierr,varid,nbdim,nid
25integer :: meanid,sdid
26integer, dimension(4)  :: id,start,sizes
27logical, save :: firstcall=.TRUE.
28integer :: l,i,j,ig0
29integer,save :: indx
30
31integer, save :: step=0
32
33! Added to work in parallel mode
34#ifdef CPP_PARA
35real px3_glop(klon_glo,llm) ! to store a 3D data set on global physics grid
36real px3_glo(iim,jjp1,llm) ! to store a global 3D data set on global lonxlat grid
37real px2_glop(klon_glo) ! to store a 2D data set on global physics grid
38real px2_glo(iim,jjp1) ! to store a 2D data set on global lonxlat grid
39real px2(ngrid)
40real px3(ngrid,llm)
41#else
42! When not running in parallel mode:
43real px3_glop(ngrid,llm) ! to store a 3D data set on global physics grid
44real px3_glo(iim,jjp1,llm) ! to store a global 3D data set on global lonxlat grid
45real px2_glop(ngrid) ! to store a 2D data set on global physics grid
46real px2_glo(iim,jjp1) ! to store a 2D data set on global lonxlat grid
47#endif
48
49! 1. Initialization (creation of stats.nc file)
50if (firstcall) then
51   firstcall=.false.
52   firstvar=trim((nom))
53   call inistats(ierr)
54endif
55
56if (firstvar==nom) then ! If we're back to the first variable, increment time counter
57      step = step + 1
58endif
59
60if (mod(step,istats).ne.0) then
61  ! if its not time to write to file, exit
62   RETURN
63endif
64
65! collect fields on a global physics grid
66#ifdef CPP_PARA
67 if (dim.eq.3) then
68  px3(1:ngrid,1:llm)=px(1:ngrid,1:llm)
69  ! Gather fieds on a "global" (without redundant longitude) array
70  call Gather(px3,px3_glop)
71!$OMP MASTER
72  if (is_mpi_root) then
73    call Grid1Dto2D_glo(px3_glop,px3_glo)
74    ! copy dx3_glo() to dx3(:) and add redundant longitude
75    dx3(1:iim,:,:)=px3_glo(1:iim,:,:)
76    dx3(iip1,:,:)=dx3(1,:,:)
77  endif
78!$OMP END MASTER
79!$OMP BARRIER
80 else ! dim.eq.2
81  ! Gather fieds on a "global" (without redundant longitude) array
82  px2(:)=px(:,1)
83  call Gather(px2,px2_glop)
84!$OMP MASTER
85          if (is_mpi_root) then
86            call Grid1Dto2D_glo(px2_glop,px2_glo)
87            ! copy px2_glo() to dx2(:) and add redundant longitude
88            dx2(1:iim,:)=px2_glo(1:iim,:)
89            dx2(iip1,:)=dx2(1,:)
90          endif
91!$OMP END MASTER
92!$OMP BARRIER
93 endif
94#else
95  if (dim.eq.3) then
96    px3_glop(:,1:llm)=px(:,1:llm)
97!  Passage variable physique -->  variable dynamique
98    DO l=1,llm
99      DO i=1,iim
100         px3_glo(i,1,l)=px(1,l)
101         px3_glo(i,jjp1,l)=px(ngrid,l)
102      ENDDO
103      DO j=2,jjm
104         ig0= 1+(j-2)*iim
105         DO i=1,iim
106            px3_glo(i,j,l)=px(ig0+i,l)
107         ENDDO
108      ENDDO
109    ENDDO
110  else ! dim.eq.2
111    px2_glop(:)=px(:,1)
112!    Passage variable physique -->  physique dynamique
113   DO i=1,iim
114     px2_glo(i,1)=px(1,1)
115     px2_glo(i,jjp1)=px(ngrid,1)
116   ENDDO
117   DO j=2,jjm
118     ig0= 1+(j-2)*iim
119     DO i=1,iim
120        px2_glo(i,j)=px(ig0+i,1)
121     ENDDO
122   ENDDO
123  endif
124#endif
125
126! 2. Write field to file
127
128if (is_master) then
129! only master needs do this
130
131ierr = NF_OPEN("stats.nc",NF_WRITE,nid)
132
133namebis=trim(nom)
134
135! test: check if that variable already exists in file
136ierr= NF_INQ_VARID(nid,namebis,meanid)
137
138if (ierr.ne.NF_NOERR) then
139  ! variable not in file, create/define it
140   if (firstvar==nom) then
141      indx=1
142      count(:)=0
143   endif
144
145
146!declaration de la variable
147
148! choix du nom des coordonnees
149   ierr= NF_INQ_DIMID(nid,"longitude",id(1))
150   ierr= NF_INQ_DIMID(nid,"latitude",id(2))
151   if (dim.eq.3) then
152      ierr= NF_INQ_DIMID(nid,"altitude",id(3))
153      ierr= NF_INQ_DIMID(nid,"Time",id(4))
154      nbdim=4
155   else if (dim==2) then
156      ierr= NF_INQ_DIMID(nid,"Time",id(3))
157      nbdim=3
158   endif
159
160   write (*,*) "====================="
161   write (*,*) "STATS: creation de ",nom
162   namebis=trim(nom)
163   call def_var(nid,namebis,titre,unite,nbdim,id,meanid,ierr)
164   if (dim.eq.3) then
165     call inivar(nid,meanid,size(px3_glop,1),dim,indx,px3_glop,ierr)
166   else ! dim.eq.2
167     call inivar(nid,meanid,size(px2_glop,1),dim,indx,px2_glop,ierr)
168   endif
169   namebis=trim(nom)//"_sd"
170   call def_var(nid,namebis,trim(titre)//" total standard deviation over the season",unite,nbdim,id,sdid,ierr)
171   if (dim.eq.3) then
172     call inivar(nid,sdid,size(px3_glop,1),dim,indx,px3_glop,ierr)
173   else ! dim.eq.2
174     call inivar(nid,sdid,size(px2_glop,1),dim,indx,px2_glop,ierr)
175   endif
176
177   ierr= NF_CLOSE(nid)
178   return
179
180else
181   ! variable found in file
182   namebis=trim(nom)//"_sd"
183   ierr= NF_INQ_VARID(nid,namebis,sdid)
184
185endif
186
187if (firstvar==nom) then
188   count(indx)=count(int(indx))+1
189   indx=indx+1
190   if (indx>istime) then
191      indx=1
192   endif
193endif
194
195if (count(indx)==0) then
196   ! very first time we write the variable to file
197   if (dim.eq.3) then
198      start=(/1,1,1,indx/)
199      sizes=(/iip1,jjp1,llm,1/)
200      mean3d(:,:,:)=0
201      sd3d(:,:,:)=0
202   else if (dim.eq.2) then
203      start=(/1,1,indx,0/)
204      sizes=(/iip1,jjp1,1,0/)
205      mean2d(:,:)=0
206      sd2d(:,:)=0
207   endif
208else
209   ! load values from file
210   if (dim.eq.3) then
211      start=(/1,1,1,indx/)
212      sizes=(/iip1,jjp1,llm,1/)
213#ifdef NC_DOUBLE
214      ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,sizes,mean3d)
215      ierr = NF_GET_VARA_DOUBLE(nid,sdid,start,sizes,sd3d)
216#else
217      ierr = NF_GET_VARA_REAL(nid,meanid,start,sizes,mean3d)
218      ierr = NF_GET_VARA_REAL(nid,sdid,start,sizes,sd3d)
219#endif
220      if (ierr.ne.NF_NOERR) then
221         write (*,*) NF_STRERROR(ierr)
222         stop ""
223      endif
224
225   else if (dim.eq.2) then
226      start=(/1,1,indx,0/)
227      sizes=(/iip1,jjp1,1,0/)
228#ifdef NC_DOUBLE
229      ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,sizes,mean2d)
230      ierr = NF_GET_VARA_DOUBLE(nid,sdid,start,sizes,sd2d)
231#else
232      ierr = NF_GET_VARA_REAL(nid,meanid,start,sizes,mean2d)
233      ierr = NF_GET_VARA_REAL(nid,sdid,start,sizes,sd2d)
234#endif
235      if (ierr.ne.NF_NOERR) then
236         write (*,*) NF_STRERROR(ierr)
237         stop ""
238      endif
239   endif
240endif ! of if (count(indx)==0)
241
242
243! 2.1. Build dx* (data on lon-lat grid, with redundant longitude)
244
245if (dim.eq.3) then
246  dx3(1:iim,:,:)=px3_glo(:,:,:)
247  dx3(iip1,:,:)=dx3(1,:,:)
248else ! dim.eq.2
249  dx2(1:iim,:)=px2_glo(:,:)
250  dx2(iip1,:)=dx2(1,:)
251endif
252
253
254! 2.2. Add current values to previously stored sums
255
256if (dim.eq.3) then
257
258   mean3d(:,:,:)=mean3d(:,:,:)+dx3(:,:,:)
259   sd3d(:,:,:)=sd3d(:,:,:)+dx3(:,:,:)**2
260
261#ifdef NC_DOUBLE
262   ierr = NF_PUT_VARA_DOUBLE(nid,meanid,start,sizes,mean3d)
263   ierr = NF_PUT_VARA_DOUBLE(nid,sdid,start,sizes,sd3d)
264#else
265   ierr = NF_PUT_VARA_REAL(nid,meanid,start,sizes,mean3d)
266   ierr = NF_PUT_VARA_REAL(nid,sdid,start,sizes,sd3d)
267#endif
268
269else if (dim.eq.2) then
270
271   mean2d(:,:)= mean2d(:,:)+dx2(:,:)
272   sd2d(:,:)=sd2d(:,:)+dx2(:,:)**2
273
274#ifdef NC_DOUBLE
275   ierr = NF_PUT_VARA_DOUBLE(nid,meanid,start,sizes,mean2d)
276   ierr = NF_PUT_VARA_DOUBLE(nid,sdid,start,sizes,sd2d)
277#else
278   ierr = NF_PUT_VARA_REAL(nid,meanid,start,sizes,mean2d)
279   ierr = NF_PUT_VARA_REAL(nid,sdid,start,sizes,sd2d)
280#endif
281
282endif ! of if (dim.eq.3) elseif (dim.eq.2)
283
284  ierr= NF_CLOSE(nid)
285endif ! of if (is_master)
286
287end subroutine wstats
288
289!======================================================
290subroutine inivar(nid,varid,ngrid,dim,indx,px,ierr)
291
292implicit none
293
294include "dimensions.h"
295!include "dimphys.h"
296include "netcdf.inc"
297
298integer, intent(in) :: nid,varid,dim,indx,ngrid
299real, dimension(ngrid,llm), intent(in) :: px
300integer, intent(out) :: ierr
301
302integer,parameter :: iip1=iim+1
303integer,parameter :: jjp1=jjm+1
304
305integer :: l,i,j,ig0
306integer, dimension(4) :: start,sizes
307real, dimension(iip1,jjp1,llm) :: dx3
308real, dimension(iip1,jjp1) :: dx2
309
310if (dim.eq.3) then
311
312   start=(/1,1,1,indx/)
313   sizes=(/iip1,jjp1,llm,1/)
314
315!  Passage variable physique -->  variable dynamique
316
317   DO l=1,llm
318      DO i=1,iip1
319         dx3(i,1,l)=px(1,l)
320         dx3(i,jjp1,l)=px(ngrid,l)
321      ENDDO
322      DO j=2,jjm
323         ig0= 1+(j-2)*iim
324         DO i=1,iim
325            dx3(i,j,l)=px(ig0+i,l)
326         ENDDO
327         dx3(iip1,j,l)=dx3(1,j,l)
328      ENDDO
329   ENDDO
330
331#ifdef NC_DOUBLE
332   ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,sizes,dx3)
333#else
334   ierr = NF_PUT_VARA_REAL(nid,varid,start,sizes,dx3)
335#endif
336
337else if (dim.eq.2) then
338
339      start=(/1,1,indx,0/)
340      sizes=(/iip1,jjp1,1,0/)
341
342!    Passage variable physique -->  physique dynamique
343
344  DO i=1,iip1
345     dx2(i,1)=px(1,1)
346     dx2(i,jjp1)=px(ngrid,1)
347  ENDDO
348  DO j=2,jjm
349     ig0= 1+(j-2)*iim
350     DO i=1,iim
351        dx2(i,j)=px(ig0+i,1)
352     ENDDO
353     dx2(iip1,j)=dx2(1,j)
354  ENDDO
355
356#ifdef NC_DOUBLE
357   ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,sizes,dx2)
358#else
359   ierr = NF_PUT_VARA_REAL(nid,varid,start,sizes,dx2)
360#endif
361
362endif
363
364end subroutine inivar
365
366!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
367
368subroutine def_var_stats(nid,name,title,units,nbdim,dimids,nvarid,ierr)
369
370! This subroutine defines variable 'name' in a (pre-existing and opened)
371! NetCDF file (known from its NetCDF ID 'nid').
372! The number of dimensions 'nbdim' of the variable, as well as the IDs of
373! corresponding dimensions must be set (in array 'dimids').
374! Upon successfull definition of the variable, 'nvarid' contains the
375! NetCDF ID of the variable.
376! The variables' attributes 'title' (Note that 'long_name' would be more
377! appropriate) and 'units' are also set.
378
379implicit none
380
381#include "netcdf.inc"
382
383integer,intent(in) :: nid ! NetCDF file ID
384character(len=*),intent(in) :: name ! the variable's name
385character(len=*),intent(in) :: title ! 'title' attribute of variable
386character(len=*),intent(in) :: units ! 'units' attribute of variable
387integer,intent(in) :: nbdim ! number of dimensions of the variable
388integer,dimension(nbdim),intent(in) :: dimids ! NetCDF IDs of the dimensions
389                                              ! the variable is defined along
390integer,intent(out) :: nvarid ! NetCDF ID of the variable
391integer,intent(out) :: ierr ! returned NetCDF staus code
392
393! 1. Switch to NetCDF define mode
394ierr=NF_REDEF(nid)
395
396! 2. Define the variable
397#ifdef NC_DOUBLE
398ierr = NF_DEF_VAR (nid,adjustl(name),NF_DOUBLE,nbdim,dimids,nvarid)
399#else
400ierr = NF_DEF_VAR (nid,adjustl(name),NF_FLOAT,nbdim,dimids,nvarid)
401#endif
402if(ierr/=NF_NOERR) then
403   write(*,*) "def_var_stats: Failed defining variable "//trim(name)
404   write(*,*) NF_STRERROR(ierr)
405   stop ""
406endif
407
408! 3. Write attributes
409ierr=NF_PUT_ATT_TEXT(nid,nvarid,"title",&
410                     len_trim(adjustl(title)),adjustl(title))
411if(ierr/=NF_NOERR) then
412   write(*,*) "def_var_stats: Failed writing title attribute for "//trim(name)
413   write(*,*) NF_STRERROR(ierr)
414   stop ""
415endif
416
417ierr=NF_PUT_ATT_TEXT(nid,nvarid,"units",&
418                     len_trim(adjustl(units)),adjustl(units))
419if(ierr/=NF_NOERR) then
420   write(*,*) "def_var_stats: Failed writing units attribute for "//trim(name)
421   write(*,*) NF_STRERROR(ierr)
422   stop ""
423endif
424
425! 4. Switch out of NetCDF define mode
426ierr = NF_ENDDEF(nid)
427
428end subroutine def_var_stats
429
Note: See TracBrowser for help on using the repository browser.