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

Last change on this file since 1292 was 1266, checked in by aslmd, 11 years ago

LMDZ.MARS
IMPORTANT CHANGE

  • Remove all reference/use of nlayermx and dimphys.h
  • Made use of automatic arrays whenever arrays are needed with dimension nlayer
  • Remove lots of obsolete reference to dimensions.h
  • Converted iono.h and param_v4.h into corresponding modules

(with embedded subroutine to allocate arrays)
(no arrays allocated if thermosphere not used)

  • Deleted param.h and put contents into module param_v4_h
  • Adapted testphys1d, newstart, etc...
  • Made DATA arrays in param_read to be initialized by subroutine

fill_data_thermos in module param_v4_h

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