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

Last change on this file was 3282, checked in by emillour, 9 months ago

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