source: trunk/LMDZ.MARS/libf/phymars/wstats_mod.F90 @ 3026

Last change on this file since 3026 was 2997, checked in by emillour, 18 months ago

Mars PCM:
Some adaptations to enable running the 1D model with XIOS.
Note that this requires compiling with "-io xios -parallel mpi"
but the model should then be run using a single core, i.e. without mpirun
EM
File size: 27.6 KB
Line 
1module wstats_mod
2
3! module containing parameters and routines to generate the "stats.nc" file
4! which will contain a mean statistical day of variables extracted at set
5! times of day every day of the simulation, and the standard deviations thereof
6
7implicit none
8
9logical,save :: callstats ! global flag to trigger generating a stats.nc
10                          ! file or not. Initialized in conf_gcm()
11!$OMP THREADPRIVATE(callstats)
12
13integer,save :: istats ! calculate stats every istats physics timestep,
14                       ! starting at first call. Initialized by inistats()
15!$OMP THREADPRIVATE(istats)
16
17integer,parameter :: istime=12 ! number of time steps per sol to store
18
19integer,save :: count(istime) ! count tab to know the variable record
20!$OMP THREADPRIVATE(count)
21
22contains
23
24subroutine wstats(ngrid,nom,titre,unite,dim,px)
25
26! Main routine to call to trigger writing a given field to the "stats.nc" file
27
28use mod_phys_lmdz_para, only : is_mpi_root, is_master, gather, klon_mpi_begin
29use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo, &
30                              nbp_lon, nbp_lat, nbp_lev, &
31                              grid_type, unstructured
32implicit none
33
34include "netcdf.inc"
35
36integer,intent(in) :: ngrid
37character (len=*),intent(in) :: nom,titre,unite
38integer,intent(in) :: dim
39real,intent(in) :: px(ngrid,nbp_lev)
40real,allocatable,save :: mean3d(:,:,:),sd3d(:,:,:),dx3(:,:,:)
41real,allocatable,save :: mean2d(:,:),sd2d(:,:),dx2(:,:)
42character (len=50) :: namebis
43character (len=50), save :: firstvar
44!$OMP THREADPRIVATE(mean3d,sd3d,dx3,mean2d,sd2d,dx2,firstvar)
45integer :: ierr,varid,nbdim,nid
46integer :: meanid,sdid
47integer, dimension(4)  :: id,start,sizes
48logical, save :: firstcall=.TRUE.
49integer,save :: indx
50integer, save :: step=0
51!$OMP THREADPRIVATE(firstcall,indx,step)
52integer :: l,i,j,ig0,n
53
54! Added to read an optional stats.def file to select outputs
55logical,save :: stats_def ! .true. if there is a stats.def file
56integer,save :: n_name_def ! number of fields to output in stats.nc
57! NB: stats_def and n_name_def do not need be threadprivate
58integer,parameter :: n_name_def_max=199 ! max number of fields to output
59character(len=120),save :: name_def(n_name_def_max)
60logical :: getout ! to trigger an early exit if variable not in output list
61
62!$OMP THREADPRIVATE(stats_def,n_name_def,name_def)
63
64! Added to work in parallel mode
65#ifdef CPP_PARA
66real px3_glop(klon_glo,nbp_lev) ! to store a 3D data set on global physics grid
67real px3_glo(nbp_lon,nbp_lat,nbp_lev) ! to store a global 3D data set on global lonxlat grid
68real px2_glop(klon_glo) ! to store a 2D data set on global physics grid
69real px2_glo(nbp_lon,nbp_lat) ! to store a 2D data set on global lonxlat grid
70real px2(ngrid)
71real px3(ngrid,nbp_lev)
72#else
73! When not running in parallel mode:
74real px3_glop(ngrid,nbp_lev) ! to store a 3D data set on global physics grid
75real px3_glo(nbp_lon,nbp_lat,nbp_lev) ! to store a global 3D data set on global lonxlat grid
76real px2_glop(ngrid) ! to store a 2D data set on global physics grid
77real px2_glo(nbp_lon,nbp_lat) ! to store a 2D data set on global lonxlat grid
78#endif
79
80! 0. Preliminary stuff
81if (callstats.eqv..false.) then
82  ! exit because creating/writing stats.nc not requested by user
83  return
84endif
85
86if (grid_type==unstructured) then
87  ! exit because non-structured grid case is not handled
88  return
89endif
90
91! 1. Initialization (creation of stats.nc file)
92if (firstcall) then
93   firstcall=.false.
94
95!$OMP MASTER
96  ! open stats.def definition file if there is one
97  open(99,file="stats.def",status='old',form='formatted',&
98       iostat=ierr)
99  if (ierr.eq.0) then
100    stats_def=.true. ! yes there is a stats.def file
101    write(*,*) "*****************"
102    write(*,*) "Reading stats.def"
103    write(*,*) "*****************"
104    do n=1,n_name_def_max
105      read(99,fmt='(a)',end=88) name_def(n)
106      write(*,*) 'Output in stats: ', trim(name_def(n))
107    enddo
10888  continue
109    ! check there is no overflow
110    if (n.ge.n_name_def_max) then
111      write(*,*) "n_name_def_max too small in wstats:",n
112      call abort_physic("wstats","n_name_def_max too small",1)
113    endif
114    n_name_def=n-1
115    close(99)
116  else
117    stats_def=.false. ! no stats.def file; output all fields sent to wstats
118  endif ! of if (ierr.eq.0)
119!$OMP END MASTER
120!$OMP BARRIER
121
122   firstvar=trim((nom))
123   call inistats(ierr)
124   if (klon_glo>1) then ! general case, 3D GCM
125     allocate(mean3d(nbp_lon+1,nbp_lat,nbp_lev))
126     allocate(sd3d(nbp_lon+1,nbp_lat,nbp_lev))
127     allocate(dx3(nbp_lon+1,nbp_lat,nbp_lev))
128     allocate(mean2d(nbp_lon+1,nbp_lat))
129     allocate(sd2d(nbp_lon+1,nbp_lat))
130     allocate(dx2(nbp_lon+1,nbp_lat))
131   else ! 1D model
132     allocate(mean3d(1,1,nbp_lev))
133     allocate(sd3d(1,1,nbp_lev))
134     allocate(dx3(1,1,nbp_lev))
135     allocate(mean2d(1,1))
136     allocate(sd2d(1,1))
137     allocate(dx2(1,1))
138   endif
139endif ! of if (firstcall)
140
141if (firstvar==nom) then ! If we're back to the first variable, increment time counter
142      step = step + 1
143endif
144
145if (mod(step,istats).ne.0) then
146  ! if its not time to write to file, exit
147   RETURN
148endif
149
150! Exit if there is a stats.def file and the variable is not in the list
151if (stats_def) then
152  getout=.true.
153  do n=1,n_name_def
154    ! look for the variable's name in the list
155    if (trim(name_def(n)).eq.nom) then
156      getout=.false.
157      ! found it, no need to scan the rest of the list exit loop
158      exit
159    endif
160  enddo
161  if (getout) then
162    ! variable not in the list so exit routine now
163    return
164  endif
165endif ! of if (stats_def)
166
167! collect fields on a global physics grid
168#ifdef CPP_PARA
169 if (dim.eq.3) then
170  px3(1:ngrid,1:nbp_lev)=px(1:ngrid,1:nbp_lev)
171  ! Gather fieds on a "global" (without redundant longitude) array
172  call Gather(px3,px3_glop)
173!$OMP MASTER
174  if (is_mpi_root) then
175    call Grid1Dto2D_glo(px3_glop,px3_glo)
176    if (klon_glo>1) then ! general case (unless in 1D)
177    ! copy dx3_glo() to dx3(:) and add redundant longitude
178    dx3(1:nbp_lon,:,:)=px3_glo(1:nbp_lon,:,:)
179    dx3(nbp_lon+1,:,:)=dx3(1,:,:)
180    endif
181  endif
182!$OMP END MASTER
183!$OMP BARRIER
184 else ! dim.eq.2
185  ! Gather fieds on a "global" (without redundant longitude) array
186  px2(:)=px(:,1)
187  call Gather(px2,px2_glop)
188!$OMP MASTER
189  if (is_mpi_root) then
190    call Grid1Dto2D_glo(px2_glop,px2_glo)
191    if (klon_glo>1) then ! general case (unless in 1D)
192      ! copy px2_glo() to dx2(:) and add redundant longitude
193      dx2(1:nbp_lon,:)=px2_glo(1:nbp_lon,:)
194      dx2(nbp_lon+1,:)=dx2(1,:)
195    endif
196  endif
197!$OMP END MASTER
198!$OMP BARRIER
199 endif
200#else
201  if (dim.eq.3) then
202    px3_glop(:,1:nbp_lev)=px(:,1:nbp_lev)
203!  Passage variable physique -->  variable dynamique
204    DO l=1,nbp_lev
205      DO i=1,nbp_lon
206         px3_glo(i,1,l)=px(1,l)
207         px3_glo(i,nbp_lat,l)=px(ngrid,l)
208      ENDDO
209      DO j=2,nbp_lat-1
210         ig0= 1+(j-2)*nbp_lon
211         DO i=1,nbp_lon
212            px3_glo(i,j,l)=px(ig0+i,l)
213         ENDDO
214      ENDDO
215    ENDDO
216  else ! dim.eq.2
217    px2_glop(:)=px(:,1)
218!    Passage variable physique -->  physique dynamique
219   DO i=1,nbp_lon
220     px2_glo(i,1)=px(1,1)
221     px2_glo(i,nbp_lat)=px(ngrid,1)
222   ENDDO
223   DO j=2,nbp_lat-1
224     ig0= 1+(j-2)*nbp_lon
225     DO i=1,nbp_lon
226        px2_glo(i,j)=px(ig0+i,1)
227     ENDDO
228   ENDDO
229  endif
230#endif
231
232! 2. Write field to file
233
234if (is_master) then
235! only master needs do this
236
237ierr = NF_OPEN("stats.nc",NF_WRITE,nid)
238
239namebis=trim(nom)
240
241! test: check if that variable already exists in file
242ierr= NF_INQ_VARID(nid,namebis,meanid)
243
244if (ierr.ne.NF_NOERR) then
245  ! variable not in file, create/define it
246   if (firstvar==nom) then
247      indx=1
248      count(:)=0
249   endif
250
251
252!declaration de la variable
253
254! choix du nom des coordonnees
255   ierr= NF_INQ_DIMID(nid,"longitude",id(1))
256   ierr= NF_INQ_DIMID(nid,"latitude",id(2))
257   if (dim.eq.3) then
258      ierr= NF_INQ_DIMID(nid,"altitude",id(3))
259      ierr= NF_INQ_DIMID(nid,"Time",id(4))
260      nbdim=4
261   else if (dim==2) then
262      ierr= NF_INQ_DIMID(nid,"Time",id(3))
263      nbdim=3
264   endif
265
266   write (*,*) "====================="
267   write (*,*) "STATS: creating ",nom
268   namebis=trim(nom)
269   call def_var_stats(nid,namebis,titre,unite,nbdim,id,meanid,ierr)
270   if (dim.eq.3) then
271     call inivar(nid,meanid,size(px3_glop,1),dim,indx,px3_glop,ierr)
272   else ! dim.eq.2
273     call inivar(nid,meanid,size(px2_glop,1),dim,indx,px2_glop,ierr)
274   endif
275   namebis=trim(nom)//"_sd"
276   call def_var_stats(nid,namebis,trim(titre)//" total standard deviation over the season",unite,nbdim,id,sdid,ierr)
277   if (dim.eq.3) then
278     call inivar(nid,sdid,size(px3_glop,1),dim,indx,px3_glop,ierr)
279   else ! dim.eq.2
280     call inivar(nid,sdid,size(px2_glop,1),dim,indx,px2_glop,ierr)
281   endif
282
283   ierr= NF_CLOSE(nid)
284   return
285
286else
287   ! variable found in file
288   namebis=trim(nom)//"_sd"
289   ierr= NF_INQ_VARID(nid,namebis,sdid)
290
291endif
292
293if (firstvar==nom) then
294   count(indx)=count(int(indx))+1
295   indx=indx+1
296   if (indx>istime) then
297      indx=1
298   endif
299endif
300
301if (count(indx)==0) then
302   ! very first time we write the variable to file
303   if (dim.eq.3) then
304      start=(/1,1,1,indx/)
305      if (klon_glo>1) then !general case
306        sizes=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
307      else
308        sizes=(/1,1,nbp_lev,1/)
309      endif
310      mean3d(:,:,:)=0
311      sd3d(:,:,:)=0
312   else if (dim.eq.2) then
313      start=(/1,1,indx,0/)
314      if (klon_glo>1) then !general case
315        sizes=(/nbp_lon+1,nbp_lat,1,0/)
316      else
317        sizes=(/1,1,1,0/)
318      endif
319      mean2d(:,:)=0
320      sd2d(:,:)=0
321   endif
322else
323   ! load values from file
324   if (dim.eq.3) then
325      start=(/1,1,1,indx/)
326      if (klon_glo>1) then !general case
327        sizes=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
328      else ! 1D model case
329        sizes=(/1,1,nbp_lev,1/)
330      endif
331#ifdef NC_DOUBLE
332      ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,sizes,mean3d)
333      ierr = NF_GET_VARA_DOUBLE(nid,sdid,start,sizes,sd3d)
334#else
335      ierr = NF_GET_VARA_REAL(nid,meanid,start,sizes,mean3d)
336      ierr = NF_GET_VARA_REAL(nid,sdid,start,sizes,sd3d)
337#endif
338      if (ierr.ne.NF_NOERR) then
339         write (*,*) "wstats error reading :",trim(nom)
340         write (*,*) NF_STRERROR(ierr)
341         call abort_physic("wstats","Failed reading "//trim(nom),1)
342      endif
343
344   else if (dim.eq.2) then
345      start=(/1,1,indx,0/)
346      if (klon_glo>1) then ! general case
347        sizes=(/nbp_lon+1,nbp_lat,1,0/)
348      else
349        sizes=(/1,1,1,0/)
350      endif
351#ifdef NC_DOUBLE
352      ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,sizes,mean2d)
353      ierr = NF_GET_VARA_DOUBLE(nid,sdid,start,sizes,sd2d)
354#else
355      ierr = NF_GET_VARA_REAL(nid,meanid,start,sizes,mean2d)
356      ierr = NF_GET_VARA_REAL(nid,sdid,start,sizes,sd2d)
357#endif
358      if (ierr.ne.NF_NOERR) then
359         write (*,*) "wstats error reading :",trim(nom)
360         write (*,*) NF_STRERROR(ierr)
361         call abort_physic("wstats","Failed reading "//trim(nom),1)
362      endif
363   endif
364endif ! of if (count(indx)==0)
365
366
367! 2.1. Build dx* (data on lon-lat grid, with redundant longitude)
368
369if (dim.eq.3) then
370  dx3(1:nbp_lon,:,:)=px3_glo(:,:,:)
371  IF (klon_glo>1) THEN ! in 3D, add redundant longitude point
372    dx3(nbp_lon+1,:,:)=dx3(1,:,:)
373  ENDIF
374else ! dim.eq.2
375  dx2(1:nbp_lon,:)=px2_glo(:,:)
376  IF (klon_glo>1) THEN ! in 3D, add redundant longitude point
377    dx2(nbp_lon+1,:)=dx2(1,:)
378  ENDIF
379endif
380
381
382! 2.2. Add current values to previously stored sums
383
384if (dim.eq.3) then
385
386   mean3d(:,:,:)=mean3d(:,:,:)+dx3(:,:,:)
387   sd3d(:,:,:)=sd3d(:,:,:)+dx3(:,:,:)**2
388
389#ifdef NC_DOUBLE
390   ierr = NF_PUT_VARA_DOUBLE(nid,meanid,start,sizes,mean3d)
391   ierr = NF_PUT_VARA_DOUBLE(nid,sdid,start,sizes,sd3d)
392#else
393   ierr = NF_PUT_VARA_REAL(nid,meanid,start,sizes,mean3d)
394   ierr = NF_PUT_VARA_REAL(nid,sdid,start,sizes,sd3d)
395#endif
396  if (ierr.ne.NF_NOERR) then
397     write (*,*) "wstats error writing :",trim(nom)
398     write (*,*) NF_STRERROR(ierr)
399     call abort_physic("wstats","Failed writing "//trim(nom),1)
400  endif
401
402else if (dim.eq.2) then
403
404   mean2d(:,:)= mean2d(:,:)+dx2(:,:)
405   sd2d(:,:)=sd2d(:,:)+dx2(:,:)**2
406
407#ifdef NC_DOUBLE
408   ierr = NF_PUT_VARA_DOUBLE(nid,meanid,start,sizes,mean2d)
409   ierr = NF_PUT_VARA_DOUBLE(nid,sdid,start,sizes,sd2d)
410#else
411   ierr = NF_PUT_VARA_REAL(nid,meanid,start,sizes,mean2d)
412   ierr = NF_PUT_VARA_REAL(nid,sdid,start,sizes,sd2d)
413#endif
414  if (ierr.ne.NF_NOERR) then
415     write (*,*) "wstats error writing :",trim(nom)
416     write(*,*) "start:",start
417     write(*,*) "sizes:",sizes
418     write(*,*) "mean2d:",mean2d
419     write(*,*) "sd2d:",sd2d
420     write (*,*) NF_STRERROR(ierr)
421     call abort_physic("wstats","Failed writing "//trim(nom),1)
422  endif
423
424endif ! of if (dim.eq.3) elseif (dim.eq.2)
425
426  ierr= NF_CLOSE(nid)
427endif ! of if (is_master)
428
429end subroutine wstats
430
431!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
432
433subroutine inistats(ierr)
434
435! routine to initialize the stats file (i.e. create file, write
436! time independent coordinates, etc.)
437
438use mod_phys_lmdz_para, only : is_master
439USE vertical_layers_mod, ONLY: ap,bp,aps,bps,preff,pseudoalt,presnivs
440USE nrtype, ONLY: pi
441USE time_phylmdz_mod, ONLY: daysec,dtphys
442USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
443USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev
444implicit none
445
446include "netcdf.inc"
447
448integer,intent(out) :: ierr
449
450integer :: nid
451integer :: l,nsteppd
452real, dimension(nbp_lev) ::  sig_s
453real,allocatable :: lon_reg_ext(:) ! extended longitudes
454integer :: idim_lat,idim_lon,idim_llm,idim_llmp1,idim_time
455real, dimension(istime) :: lt
456integer :: nvarid
457
458
459IF (nbp_lon*nbp_lat==1) THEN
460  ! 1D model
461  ALLOCATE(lon_reg_ext(1))
462ELSE
463  ! 3D model
464  ALLOCATE(lon_reg_ext(nbp_lon+1))
465ENDIF
466
467write (*,*)
468write (*,*) '                        || STATS ||'
469write (*,*)
470write (*,*) 'daysec',daysec
471write (*,*) 'dtphys',dtphys
472nsteppd=nint(daysec/dtphys)
473write (*,*) 'nsteppd=',nsteppd
474
475if (abs(float(nsteppd)-daysec/dtphys).gt.1.e-8*daysec) then
476  call abort_physic("inistats",'1 sol .ne. n physics steps',1)
477endif
478
479if(mod(nsteppd,istime).ne.0) then
480  call abort_physic("inistats",'1 sol .ne. n*istime physics steps',1)
481endif
482
483istats=nsteppd/istime
484write (*,*) 'istats=',istats
485write (*,*) 'Storing ',istime,'times per day'
486write (*,*) 'thus every ',istats,'physical timestep '
487write (*,*)
488
489do l= 1, nbp_lev
490  sig_s(l)=((ap(l)+ap(l+1))/preff+bp(l)+bp(l+1))/2.
491  pseudoalt(l)=-10.*log(presnivs(l)/preff)   
492enddo
493
494lon_reg_ext(1:nbp_lon)=lon_reg(1:nbp_lon)
495IF (nbp_lon*nbp_lat/=1) THEN
496  ! In 3D, add extra redundant point (180 degrees,
497  ! since lon_reg starts at -180)
498  lon_reg_ext(nbp_lon+1)=-lon_reg_ext(1)
499ENDIF
500
501if (is_master) then
502  ! only the master needs do this
503  ierr = NF_CREATE("stats.nc",IOR(NF_CLOBBER,NF_64BIT_OFFSET),nid)
504  if (ierr.ne.NF_NOERR) then
505    write (*,*) NF_STRERROR(ierr)
506    call abort_physic("inistats",'failed creating stats.nc',1)
507  endif
508
509  ierr = NF_DEF_DIM (nid, "latitude", nbp_lat, idim_lat)
510  IF (nbp_lon*nbp_lat==1) THEN
511    ierr = NF_DEF_DIM (nid, "longitude", 1, idim_lon)
512  ELSE
513    ierr = NF_DEF_DIM (nid, "longitude", nbp_lon+1, idim_lon)
514  ENDIF
515  ierr = NF_DEF_DIM (nid, "altitude", nbp_lev, idim_llm)
516  ierr = NF_DEF_DIM (nid, "llmp1", nbp_lev+1, idim_llmp1)
517  ierr = NF_DEF_DIM (nid, "Time", NF_UNLIMITED, idim_time)
518
519  ierr = NF_ENDDEF(nid)
520
521  call def_var_stats(nid,"Time","Time", &
522                     "days since 0000-00-0 00:00:00",1, &
523                     [idim_time],nvarid,ierr)
524  ! Time is initialised later by mkstats subroutine
525
526  call def_var_stats(nid,"latitude","latitude", &
527                     "degrees_north",1,[idim_lat],nvarid,ierr)
528#ifdef NC_DOUBLE
529  ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lat_reg/pi*180)
530#else
531  ierr = NF_PUT_VAR_REAL (nid,nvarid,lat_reg/pi*180)
532#endif
533
534  call def_var_stats(nid,"longitude","East longitude", &
535                     "degrees_east",1,[idim_lon],nvarid,ierr)
536#ifdef NC_DOUBLE
537  ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lon_reg_ext/pi*180)
538#else
539  ierr = NF_PUT_VAR_REAL (nid,nvarid,lon_reg_ext/pi*180)
540#endif
541
542  ! Niveaux verticaux, aps et bps
543  ierr = NF_REDEF (nid)
544#ifdef NC_DOUBLE
545  ierr = NF_DEF_VAR (nid,"altitude", NF_DOUBLE, 1,[idim_llm],nvarid)
546#else
547  ierr = NF_DEF_VAR (nid,"altitude", NF_FLOAT, 1,[idim_llm],nvarid)
548#endif
549  ierr = NF_PUT_ATT_TEXT (nid,nvarid,"long_name",8,"altitude")
550  ierr = NF_PUT_ATT_TEXT (nid,nvarid,'units',2,"km")
551  ierr = NF_PUT_ATT_TEXT (nid,nvarid,'positive',2,"up")
552  ierr = NF_ENDDEF(nid)
553#ifdef NC_DOUBLE
554  ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,pseudoalt)
555#else
556  ierr = NF_PUT_VAR_REAL (nid,nvarid,pseudoalt)
557#endif
558  call def_var_stats(nid,"aps","hybrid pressure at midlayers", &
559                     " ",1,[idim_llm],nvarid,ierr)
560#ifdef NC_DOUBLE
561  ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,aps)
562#else
563  ierr = NF_PUT_VAR_REAL (nid,nvarid,aps)
564#endif
565
566  call def_var_stats(nid,"bps","hybrid sigma at midlayers", &
567                     " ",1,[idim_llm],nvarid,ierr)
568#ifdef NC_DOUBLE
569  ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,bps)
570#else
571  ierr = NF_PUT_VAR_REAL (nid,nvarid,bps)
572#endif
573
574  ierr=NF_CLOSE(nid)
575
576endif ! of if (is_master)
577
578end subroutine inistats
579
580!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
581
582subroutine inivar(nid,varid,ngrid,dim,indx,px,ierr)
583
584! routine to initialize writing a variable in the stats file
585
586use mod_grid_phy_lmdz, only : nbp_lon, nbp_lat, nbp_lev, klon_glo
587
588implicit none
589
590include "netcdf.inc"
591
592integer, intent(in) :: nid,varid,dim,indx,ngrid
593real, dimension(ngrid,nbp_lev), intent(in) :: px
594integer, intent(out) :: ierr
595
596integer :: l,i,j,ig0
597integer, dimension(4) :: start,sizes
598real, dimension(nbp_lon+1,nbp_lat,nbp_lev) :: dx3
599real, dimension(nbp_lon+1,nbp_lat) :: dx2
600real :: dx3_1d(nbp_lev) ! for 1D outputs
601real :: dx2_1d ! for 1D outputs
602
603if (dim.eq.3) then
604
605   start=(/1,1,1,indx/)
606   if (klon_glo>1) then ! general 3D case
607     sizes=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
608   else
609     sizes=(/1,1,nbp_lev,1/)
610   endif
611
612!  Passage variable physique -->  variable dynamique
613
614   if (klon_glo>1) then ! general case
615    DO l=1,nbp_lev
616      DO i=1,nbp_lon+1
617         dx3(i,1,l)=px(1,l)
618         dx3(i,nbp_lat,l)=px(ngrid,l)
619      ENDDO
620      DO j=2,nbp_lat-1
621         ig0= 1+(j-2)*nbp_lon
622         DO i=1,nbp_lon
623            dx3(i,j,l)=px(ig0+i,l)
624         ENDDO
625         dx3(nbp_lon+1,j,l)=dx3(1,j,l)
626      ENDDO
627    ENDDO
628   else ! 1D model case
629     dx3_1d(1:nbp_lev)=px(1,1:nbp_lev)
630   endif
631
632#ifdef NC_DOUBLE
633   if (klon_glo>1) then
634     ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,sizes,dx3)
635   else
636     ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,sizes,dx3_1d)
637   endif
638#else
639   if (klon_glo>1) then
640     ierr = NF_PUT_VARA_REAL(nid,varid,start,sizes,dx3)
641   else
642     ierr = NF_PUT_VARA_REAL(nid,varid,start,sizes,dx3_1d)
643   endif
644#endif
645  if (ierr.ne.NF_NOERR) then
646     write (*,*) "inivar error writing variable"
647     write (*,*) NF_STRERROR(ierr)
648     call abort_physic("inivar","error writing variable",1)
649  endif
650
651else if (dim.eq.2) then
652
653      start=(/1,1,indx,0/)
654      if (klon_glo>1) then ! general 3D case
655        sizes=(/nbp_lon+1,nbp_lat,1,0/)
656      else
657        sizes=(/1,1,1,0/)
658      endif
659
660!    Passage variable physique -->  physique dynamique
661
662  if (klon_glo>1) then ! general case
663  DO i=1,nbp_lon+1
664     dx2(i,1)=px(1,1)
665     dx2(i,nbp_lat)=px(ngrid,1)
666  ENDDO
667  DO j=2,nbp_lat-1
668     ig0= 1+(j-2)*nbp_lon
669     DO i=1,nbp_lon
670        dx2(i,j)=px(ig0+i,1)
671     ENDDO
672     dx2(nbp_lon+1,j)=dx2(1,j)
673  ENDDO
674  else ! 1D model case
675    dx2_1d=px(1,1)
676  endif
677 
678#ifdef NC_DOUBLE
679   if (klon_glo>1) then
680     ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,sizes,dx2)
681   else
682     ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,sizes,dx2_1d)
683   endif
684#else
685   if (klon_glo>1) then
686     ierr = NF_PUT_VARA_REAL(nid,varid,start,sizes,dx2)
687   else
688     ierr = NF_PUT_VARA_REAL(nid,varid,start,sizes,dx2_1d)
689   endif
690#endif
691  if (ierr.ne.NF_NOERR) then
692     write (*,*) "inivar error writing variable"
693     write (*,*) NF_STRERROR(ierr)
694     call abort_physic("inivar","error writing variable",1)
695  endif
696
697endif
698
699end subroutine inivar
700
701!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
702
703subroutine def_var_stats(nid,name,title,units,nbdim,dimids,nvarid,ierr)
704
705! This subroutine defines variable 'name' in a (pre-existing and opened)
706! NetCDF file (known from its NetCDF ID 'nid').
707! The number of dimensions 'nbdim' of the variable, as well as the IDs of
708! corresponding dimensions must be set (in array 'dimids').
709! Upon successfull definition of the variable, 'nvarid' contains the
710! NetCDF ID of the variable.
711! The variables' attributes 'title' (Note that 'long_name' would be more
712! appropriate) and 'units' are also set.
713
714implicit none
715
716include "netcdf.inc"
717
718integer,intent(in) :: nid ! NetCDF file ID
719character(len=*),intent(in) :: name ! the variable's name
720character(len=*),intent(in) :: title ! 'title' attribute of variable
721character(len=*),intent(in) :: units ! 'units' attribute of variable
722integer,intent(in) :: nbdim ! number of dimensions of the variable
723integer,dimension(nbdim),intent(in) :: dimids ! NetCDF IDs of the dimensions
724                                              ! the variable is defined along
725integer,intent(out) :: nvarid ! NetCDF ID of the variable
726integer,intent(out) :: ierr ! returned NetCDF staus code
727
728! 1. Switch to NetCDF define mode
729ierr=NF_REDEF(nid)
730
731! 2. Define the variable
732#ifdef NC_DOUBLE
733ierr = NF_DEF_VAR (nid,adjustl(name),NF_DOUBLE,nbdim,dimids,nvarid)
734#else
735ierr = NF_DEF_VAR (nid,adjustl(name),NF_FLOAT,nbdim,dimids,nvarid)
736#endif
737if(ierr/=NF_NOERR) then
738   write(*,*) "def_var_stats: Failed defining variable "//trim(name)
739   write(*,*) NF_STRERROR(ierr)
740   call abort_physic("def_var_stats","Failed defining "//trim(name),1)
741endif
742
743! 3. Write attributes
744ierr=NF_PUT_ATT_TEXT(nid,nvarid,"title",&
745                     len_trim(adjustl(title)),adjustl(title))
746if(ierr/=NF_NOERR) then
747   write(*,*) "def_var_stats: Failed writing title attribute for "//trim(name)
748   write(*,*) NF_STRERROR(ierr)
749   call abort_physic("def_var_stats","Failed writing title for "//trim(name),1)
750endif
751
752ierr=NF_PUT_ATT_TEXT(nid,nvarid,"units",&
753                     len_trim(adjustl(units)),adjustl(units))
754if(ierr/=NF_NOERR) then
755   write(*,*) "def_var_stats: Failed writing units attribute for "//trim(name)
756   write(*,*) NF_STRERROR(ierr)
757   call abort_physic("def_var_stats","Failed writing units for "//trim(name),1)
758endif
759
760! 4. Switch out of NetCDF define mode
761ierr = NF_ENDDEF(nid)
762
763end subroutine def_var_stats
764
765!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
766
767subroutine mkstats(ierr)
768
769!  This routine finalizes tha stats.nc file from sums and sums of squares
770!  to means and standard deviations. The data file is overwritten in place.
771
772use mod_phys_lmdz_para, only : is_master
773use mod_grid_phy_lmdz, only : nbp_lon, nbp_lat, nbp_lev, klon_glo
774
775implicit none
776
777include "netcdf.inc"
778
779integer,intent(out) :: ierr ! netCDF return error code
780integer :: nid,nbvar,i,ndims,lt,nvarid
781integer, dimension(4) :: id,varid,start,size
782integer, dimension(5) :: dimids
783character (len=50) :: name,nameout,units,title
784real,allocatable :: sum3d(:,:,:),square3d(:,:,:),mean3d(:,:,:),sd3d(:,:,:)
785real,allocatable :: sum2d(:,:),square2d(:,:),mean2d(:,:),sd2d(:,:)
786real, dimension(istime) :: time
787real, dimension(nbp_lat) :: lat
788real,allocatable :: lon(:)
789real, dimension(nbp_lev) :: alt
790logical :: lcopy=.true.
791!integer :: latid,lonid,altid,timeid
792integer :: meanid,sdid
793!integer, dimension(4) :: dimout
794
795if (is_master) then
796! only the master needs do all this
797
798! Incrementation of count for the last step, which is not done in wstats
799count(istime)=count(istime)+1
800
801if (klon_glo>1) then
802  allocate(lon(nbp_lon+1))
803  allocate(sum3d(nbp_lon+1,nbp_lat,nbp_lev))
804  allocate(square3d(nbp_lon+1,nbp_lat,nbp_lev))
805  allocate(mean3d(nbp_lon+1,nbp_lat,nbp_lev))
806  allocate(sd3d(nbp_lon+1,nbp_lat,nbp_lev))
807  allocate(sum2d(nbp_lon+1,nbp_lat))
808  allocate(square2d(nbp_lon+1,nbp_lat))
809  allocate(mean2d(nbp_lon+1,nbp_lat))
810  allocate(sd2d(nbp_lon+1,nbp_lat))
811else ! 1D model case
812  allocate(lon(1))
813  allocate(sum3d(1,1,nbp_lev))
814  allocate(square3d(1,1,nbp_lev))
815  allocate(mean3d(1,1,nbp_lev))
816  allocate(sd3d(1,1,nbp_lev))
817  allocate(sum2d(1,1))
818  allocate(square2d(1,1))
819  allocate(mean2d(1,1))
820  allocate(sd2d(1,1))
821endif
822
823ierr = NF_OPEN("stats.nc",NF_WRITE,nid)
824
825! We catch the id of dimensions of the stats file
826
827ierr= NF_INQ_DIMID(nid,"latitude",id(1))
828ierr= NF_INQ_DIMID(nid,"longitude",id(2))
829ierr= NF_INQ_DIMID(nid,"altitude",id(3))
830ierr= NF_INQ_DIMID(nid,"Time",id(4))
831
832ierr= NF_INQ_VARID(nid,"latitude",varid(1))
833ierr= NF_INQ_VARID(nid,"longitude",varid(2))
834ierr= NF_INQ_VARID(nid,"altitude",varid(3))
835ierr= NF_INQ_VARID(nid,"Time",varid(4))
836
837! Time initialisation
838
839do i=1,istime
840   time(i)=i*24./istime
841#ifdef NC_DOUBLE
842   ierr= NF_PUT_VARA_DOUBLE(nid,varid(4),[i],[1],time(i))
843#else
844   ierr= NF_PUT_VARA_REAL(nid,varid(4),[i],[1],time(i))
845#endif
846enddo
847
848! We catche the values of the variables
849
850#ifdef NC_DOUBLE
851         ierr = NF_GET_VAR_DOUBLE(nid,varid(1),lat)
852         ierr = NF_GET_VAR_DOUBLE(nid,varid(2),lon)
853         ierr = NF_GET_VAR_DOUBLE(nid,varid(3),alt)
854#else
855         ierr = NF_GET_VAR_REAL(nid,varid(1),lat)
856         ierr = NF_GET_VAR_REAL(nid,varid(2),lon)
857         ierr = NF_GET_VAR_REAL(nid,varid(3),alt)
858#endif
859
860! We catch the number of variables in the stats file
861ierr = NF_INQ_NVARS(nid,nbvar)
862
863! to catche the "real" number of variables (without the "additionnal variables")
864nbvar=(nbvar-4)/2
865
866do i=1,nbvar
867   nvarid=(i-1)*2+5
868
869   ! What's the variable's name?
870   ierr=NF_INQ_VARNAME(nid,nvarid,name)
871   write(*,*) "OK variable ",name
872   ! Its units?
873   units=" "
874   ierr=NF_GET_ATT_TEXT(nid,nvarid,"units",units)
875   ! Its title?
876   title=" "
877   ierr=NF_GET_ATT_TEXT(nid,nvarid,"title",title)
878   ! Its number of dimensions?   
879   ierr=NF_INQ_VARNDIMS(nid,nvarid,ndims)
880   ! Its values?
881
882   if(ndims==4) then ! lat, lon, alt & time
883
884!      dimout(1)=lonid
885!      dimout(2)=latid
886!      dimout(3)=altid
887!      dimout(4)=timeid
888
889      if (klon_glo>1) then ! general case
890        size=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
891      else ! 1D model
892        size=(/1,1,nbp_lev,1/)
893      endif
894      do lt=1,istime
895         start=(/1,1,1,lt/)
896         ! Extraction of the "source" variables
897#ifdef NC_DOUBLE
898         ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,size,sum3d)
899         ierr = NF_GET_VARA_DOUBLE(nid,nvarid+1,start,size,square3d)
900#else
901         ierr = NF_GET_VARA_REAL(nid,nvarid,start,size,sum3d)
902         ierr = NF_GET_VARA_REAL(nid,nvarid+1,start,size,square3d)
903#endif
904         ! Calculation of these nvariables
905         mean3d=sum3d/count(lt)
906         sd3d=sqrt(max(0.,square3d/count(lt)-mean3d**2))
907         ! Writing of the variables
908#ifdef NC_DOUBLE
909         ierr = NF_PUT_VARA_DOUBLE(nid,nvarid,start,size,mean3d)
910         ierr = NF_PUT_VARA_DOUBLE(nid,nvarid+1,start,size,sd3d)
911#else
912         ierr = NF_PUT_VARA_REAL(nid,nvarid,start,size,mean3d)
913         ierr = NF_PUT_VARA_REAL(nid,nvarid+1,start,size,sd3d)
914#endif
915      enddo
916
917    else if (ndims.eq.3) then
918
919!      dimout(1)=lonid
920!      dimout(2)=latid
921!      dimout(3)=timeid
922
923      if (klon_glo>1) then ! general case
924        size=(/nbp_lon+1,nbp_lat,1,0/)
925      else
926        size=(/1,1,1,0/)
927      endif
928      do lt=1,istime
929         start=(/1,1,lt,0/)
930         ! Extraction of the "source" variables
931#ifdef NC_DOUBLE
932         ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,size,sum2d)
933         ierr = NF_GET_VARA_DOUBLE(nid,nvarid+1,start,size,square2d)
934#else
935         ierr = NF_GET_VARA_REAL(nid,nvarid,start,size,sum2d)
936         ierr = NF_GET_VARA_REAL(nid,nvarid+1,start,size,square2d)
937#endif
938         ! Calculation of these variables
939         mean2d=sum2d/count(lt)
940         sd2d=sqrt(max(0.,square2d/count(lt)-mean2d**2))
941         ! Writing of the variables
942#ifdef NC_DOUBLE
943         ierr = NF_PUT_VARA_DOUBLE(nid,nvarid,start,size,mean2d)
944         ierr = NF_PUT_VARA_DOUBLE(nid,nvarid+1,start,size,sd2d)
945#else
946         ierr = NF_PUT_VARA_REAL(nid,nvarid,start,size,mean2d)
947         ierr = NF_PUT_VARA_REAL(nid,nvarid+1,start,size,sd2d)
948#endif
949      enddo
950
951    endif
952enddo
953
954ierr= NF_CLOSE(nid)
955
956endif ! of if (is_master)
957
958end subroutine mkstats
959
960!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
961
962end module wstats_mod
Note: See TracBrowser for help on using the repository browser.