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

Last change on this file since 1448 was 1422, checked in by milmd, 10 years ago

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