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

Last change on this file since 1769 was 1689, checked in by emillour, 8 years ago

Mars GCM:
Fixing a big bug (dating from revision 1528) in wstats.
EM

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