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

Last change on this file since 1330 was 1315, checked in by milmd, 10 years ago

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