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

Last change on this file since 2800 was 2616, checked in by romain.vande, 3 years ago

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