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

Last change on this file since 2613 was 2573, checked in by emillour, 3 years ago

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