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

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