source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/writediagfi.F @ 1057

Last change on this file since 1057 was 313, checked in by ymipsl, 10 years ago
  • implement splitting of XIOS file for lmdz physics
  • Termination is done properly in parallel by calling MPI_ABORT instead of abort or stop

YM

File size: 17.7 KB
Line 
1      subroutine writediagfi(ngrid,nom,titre,unite,dim,px)
2
3!  Ecriture de variables diagnostiques au choix dans la physique
4!  dans un fichier NetCDF nomme  'diagfi'. Ces variables peuvent etre
5!  3d (ex : temperature), 2d (ex : temperature de surface), ou
6!  0d (pour un scalaire qui ne depend que du temps : ex : la longitude
7!  solaire)
8!  (ou encore 1d, dans le cas de testphys1d, pour sortir une colonne)
9!  La periode d'ecriture est donnee par
10!  "ecritphy " regle dans le fichier de controle de run :  run.def
11!
12!    writediagfi peut etre appele de n'importe quelle subroutine
13!    de la physique, plusieurs fois. L'initialisation et la creation du
14!    fichier se fait au tout premier appel.
15!
16! WARNING : les variables dynamique (u,v,t,q,ps)
17!  sauvees par writediagfi avec une
18! date donnee sont legerement differentes que dans le fichier histoire car
19! on ne leur a pas encore ajoute de la dissipation et de la physique !!!
20! IL est  RECOMMANDE d'ajouter les tendance physique a ces variables
21! avant l'ecriture dans diagfi (cf. physiq.F)
22
23! Modifs: Aug.2010 Ehouarn: enforce outputs to be real*4
24!         Oct 2011 Francois: enable having a 'diagfi.def' file to select
25!                            at runtime, which variables to put in file
26!
27!  parametres (input) :
28!  ----------
29!      ngrid : nombres de point ou est calcule la physique
30!                (ngrid = 2+(jjm-1)*iim - 1/jjm)
31!                 (= nlon ou klon dans la physique terrestre)
32!     
33!      unit : unite logique du fichier de sortie (toujours la meme)
34!      nom  : nom de la variable a sortir (chaine de caracteres)
35!      titre: titre de la variable (chaine de caracteres)
36!      unite : unite de la variable (chaine de caracteres)
37!      px : variable a sortir (real 0, 1, 2, ou 3d)
38!      dim : dimension de px : 0, 1, 2, ou 3 dimensions
39!
40!=================================================================
41      use surfdat_h, only: phisfi
42      use control_mod, only: ecritphy, day_step, iphysiq
43      USE mod_phys_lmdz_para, only : is_parallel, is_mpi_master,
44     &                               is_master, gather
45      USE mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
46      implicit none
47
48! Commons
49      include "dimensions.h"
50      include "paramet.h"
51      include "comvert.h"
52      include "comgeom.h"
53      include "netcdf.inc"
54      include "temps.h"
55
56! Arguments on input:
57      integer,intent(in) :: ngrid
58      character (len=*),intent(in) :: nom,titre,unite
59      integer,intent(in) :: dim
60      real,intent(in) :: px(ngrid,llm)
61
62! Local variables:
63
64      real*4 dx3(iip1,jjp1,llm) ! to store a 3D data set
65      real*4 dx2(iip1,jjp1)     ! to store a 2D (surface) data set
66      real*4 dx1(llm)           ! to store a 1D (column) data set
67      real*4 dx0
68
69      real*4,save :: date
70!$OMP THREADPRIVATE(date)
71
72      REAL phis(ip1jmp1)
73
74      integer irythme
75      integer ierr,ierr2
76      integer iq
77      integer i,j,l,zmax , ig0
78
79      integer,save :: zitau=0
80      character(len=20),save :: firstnom='1234567890'
81!$OMP THREADPRIVATE(zitau,firstnom)
82
83! Ajouts
84      integer, save :: ntime=0
85!$OMP THREADPRIVATE(ntime)
86      integer :: idim,varid
87      integer :: nid
88      character(len=*),parameter :: fichnom="diagfi.nc"
89      integer, dimension(4) :: id
90      integer, dimension(4) :: edges,corner
91
92! Added to use diagfi.def to select output variable
93      logical,save :: diagfi_def
94      logical :: getout
95      integer,save :: n_nom_def
96      integer :: n
97      integer,parameter :: n_nom_def_max=199
98      character(len=120),save :: nom_def(n_nom_def_max)
99      logical,save :: firstcall=.true.
100!$OMP THREADPRIVATE(firstcall)  !diagfi_def,n_nom_def,nom_def read in diagfi.def
101     
102#ifndef MESOSCALE
103
104#ifdef CPP_PARA
105! Added to work in parallel mode
106      real dx3_glop(klon_glo,llm)
107      real dx3_glo(iim,jjp1,llm) ! to store a global 3D data set
108      real dx2_glop(klon_glo)
109      real dx2_glo(iim,jjp1)     ! to store a global 2D (surface) data set
110      real px2(ngrid)
111!      real dx1_glo(llm)          ! to store a 1D (column) data set
112!      real dx0_glo
113      real phisfi_glo(klon_glo) ! surface geopotential on global physics grid
114#else
115      real phisfi_glo(ngrid) ! surface geopotential on global physics grid
116#endif
117     
118      RETURN
119!***************************************************************
120!Sortie des variables au rythme voulu
121
122      irythme = ecritphy ! sortie au rythme de ecritphy
123!     irythme = iconser  ! sortie au rythme des variables de controle
124!     irythme = iphysiq  ! sortie a tous les pas physique
125!     irythme = iecri*day_step ! sortie au rythme des fichiers histoires
126!     irythme = periodav*day_step ! sortie au rythme des fichiers histmoy
127
128!***************************************************************
129
130! At very first call, check if there is a "diagfi.def" to use and read it
131! -----------------------------------------------------------------------
132      IF (firstcall) THEN
133         firstcall=.false.
134
135!$OMP MASTER
136  !      Open diagfi.def definition file if there is one:
137         open(99,file="diagfi.def",status='old',form='formatted',
138     s   iostat=ierr2)
139
140         if(ierr2.eq.0) then
141            diagfi_def=.true.
142            write(*,*) "******************"
143            write(*,*) "Reading diagfi.def"
144            write(*,*) "******************"
145            do n=1,n_nom_def_max
146              read(99,fmt='(a)',end=88) nom_def(n)
147              write(*,*) 'Output in diagfi: ', trim(nom_def(n))
148            end do
149 88         continue
150            if (n.ge.n_nom_def_max) then
151               write(*,*)"n_nom_def_max too small in writediagfi.F:",n
152               CALL abort_physiq
153            end if
154            n_nom_def=n-1
155            close(99)
156         else
157            diagfi_def=.false.
158         endif
159!$OMP END MASTER
160!$OMP BARRIER
161      END IF ! of IF (firstcall)
162
163! Get out of write_diagfi if there is diagfi.def AND variable not listed
164!  ---------------------------------------------------------------------
165      if (diagfi_def) then
166          getout=.true.
167          do n=1,n_nom_def
168             if(trim(nom_def(n)).eq.nom) getout=.false.
169          end do
170          if (getout) return
171      end if
172
173! Initialisation of 'firstnom' and create/open the "diagfi.nc" NetCDF file
174! ------------------------------------------------------------------------
175! (at very first call to the subroutine, in accordance with diagfi.def)
176
177      if (firstnom.eq.'1234567890') then ! .true. for the very first valid
178      !   call to this subroutine; now set 'firstnom'
179         firstnom = nom
180         ! just to be sure, check that firstnom is large enough to hold nom
181         if (len_trim(firstnom).lt.len_trim(nom)) then
182           write(*,*) "writediagfi: Error !!!"
183           write(*,*) "   firstnom string not long enough!!"
184           write(*,*) "   increase its size to at least ",len_trim(nom)
185           CALL abort_physiq
186         endif
187         
188         zitau = -1 ! initialize zitau
189
190#ifdef CPP_PARA
191          ! Gather phisfi() geopotential on physics grid
192          call Gather(phisfi,phisfi_glo)
193#else
194         phisfi_glo(:)=phisfi(:)
195#endif
196
197         !! parallel: we cannot use the usual writediagfi method
198!!         call iophys_ini
199         if (is_master) then
200         ! only the master is required to do this
201
202         ! Create the NetCDF file
203         ierr = NF_CREATE(fichnom, NF_CLOBBER, nid)
204         ! Define the 'Time' dimension
205         ierr = nf_def_dim(nid,"Time",NF_UNLIMITED,idim)
206         ! Define the 'Time' variable
207!#ifdef NC_DOUBLE
208!         ierr = NF_DEF_VAR (nid, "Time", NF_DOUBLE, 1, idim,varid)
209!#else
210         ierr = NF_DEF_VAR (nid, "Time", NF_FLOAT, 1, idim,varid)
211!#endif
212         ! Add a long_name attribute
213         ierr = NF_PUT_ATT_TEXT (nid, varid, "long_name",
214     .          4,"Time")
215         ! Add a units attribute
216         ierr = NF_PUT_ATT_TEXT(nid, varid,'units',29,
217     .          "days since 0000-00-0 00:00:00")
218         ! Switch out of NetCDF Define mode
219         ierr = NF_ENDDEF(nid)
220
221         ! write "header" of file (longitudes, latitudes, geopotential, ...)
222         call gr_fi_dyn(1,size(phisfi_glo),iip1,jjp1,phisfi_glo,phis)
223         call iniwrite(nid,day_ini,phis)
224
225         endif ! of if (is_master)
226
227      else
228
229         if (is_master) then
230           ! only the master is required to do this
231
232           ! Open the NetCDF file
233           ierr = NF_OPEN(fichnom,NF_WRITE,nid)
234         endif ! of if (is_master)
235
236      endif ! if (firstnom.eq.'1234567890')
237
238      if (ngrid.eq.1) then
239        ! in testphys1d, for the 1d version of the GCM, iphysiq and irythme
240        ! are undefined; so set them to 1
241        iphysiq=1 
242        irythme=1
243        ! NB:
244      endif
245
246! Increment time index 'zitau' if it is the "fist call" (at given time level)
247! to writediagfi
248!------------------------------------------------------------------------
249      if (nom.eq.firstnom) then
250          zitau = zitau + iphysiq
251      end if
252
253!--------------------------------------------------------
254! Write the variables to output file if it's time to do so
255!--------------------------------------------------------
256
257      if ( MOD(zitau+1,irythme) .eq.0.) then
258
259! Compute/write/extend 'Time' coordinate (date given in days)
260! (done every "first call" (at given time level) to writediagfi)
261! Note: date is incremented as 1 step ahead of physics time
262!--------------------------------------------------------
263
264        if (is_master) then
265           ! only the master is required to do this
266        if (nom.eq.firstnom) then
267        ! We have identified a "first call" (at given date)
268           ntime=ntime+1 ! increment # of stored time steps
269           ! compute corresponding date (in days and fractions thereof)
270           date= float (zitau +1)/float (day_step)
271           ! Get NetCDF ID of 'Time' variable
272           ierr= NF_INQ_VARID(nid,"Time",varid)
273           ! Write (append) the new date to the 'Time' array
274!#ifdef NC_DOUBLE
275!           ierr= NF_PUT_VARA_DOUBLE(nid,varid,ntime,1,date)
276!#else
277           ierr= NF_PUT_VARA_REAL(nid,varid,ntime,1,date)
278!#endif
279           if (ierr.ne.NF_NOERR) then
280              write(*,*) "***** PUT_VAR matter in writediagfi_nc"
281              write(*,*) "***** with time"
282              write(*,*) 'ierr=', ierr   
283c             call abort
284           endif
285
286           write(6,*)'WRITEDIAGFI: date= ', date
287        end if ! of if (nom.eq.firstnom)
288
289        endif ! of if (is_master)
290
291!Case of a 3D variable
292!---------------------
293        if (dim.eq.3) then
294
295#ifdef CPP_PARA
296          ! Gather field on a "global" (without redundant longitude) array
297          call Gather(px,dx3_glop)
298!$OMP MASTER
299          if (is_mpi_master) then
300            call Grid1Dto2D_glo(dx3_glop,dx3_glo)
301            ! copy dx3_glo() to dx3(:) and add redundant longitude
302            dx3(1:iim,:,:)=dx3_glo(1:iim,:,:)
303            dx3(iip1,:,:)=dx3(1,:,:)
304          endif
305!$OMP END MASTER
306!$OMP BARRIER
307#else
308!         Passage variable physique -->  variable dynamique
309!         recast (copy) variable from physics grid to dynamics grid
310           DO l=1,llm
311             DO i=1,iip1
312                dx3(i,1,l)=px(1,l)
313                dx3(i,jjp1,l)=px(ngrid,l)
314             ENDDO
315             DO j=2,jjm
316                ig0= 1+(j-2)*iim
317                DO i=1,iim
318                   dx3(i,j,l)=px(ig0+i,l)
319                ENDDO
320                dx3(iip1,j,l)=dx3(1,j,l)
321             ENDDO
322           ENDDO
323#endif
324!         Ecriture du champs
325
326          if (is_master) then
327           ! only the master writes to output
328! name of the variable
329           ierr= NF_INQ_VARID(nid,nom,varid)
330           if (ierr /= NF_NOERR) then
331! corresponding dimensions
332              ierr= NF_INQ_DIMID(nid,"longitude",id(1))
333              ierr= NF_INQ_DIMID(nid,"latitude",id(2))
334              ierr= NF_INQ_DIMID(nid,"altitude",id(3))
335              ierr= NF_INQ_DIMID(nid,"Time",id(4))
336
337! Create the variable if it doesn't exist yet
338
339              write (*,*) "=========================="
340              write (*,*) "DIAGFI: creating variable ",nom
341              call def_var(nid,nom,titre,unite,4,id,varid,ierr)
342
343           endif
344
345           corner(1)=1
346           corner(2)=1
347           corner(3)=1
348           corner(4)=ntime
349
350           edges(1)=iip1
351           edges(2)=jjp1
352           edges(3)=llm
353           edges(4)=1
354!#ifdef NC_DOUBLE
355!           ierr= NF_PUT_VARA_DOUBLE(nid,varid,corner,edges,dx3)
356!#else
357!           write(*,*)"test:  nid=",nid," varid=",varid
358!           write(*,*)"       corner()=",corner
359!           write(*,*)"       edges()=",edges
360!           write(*,*)"       dx3()=",dx3
361           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx3)
362!#endif
363
364           if (ierr.ne.NF_NOERR) then
365              write(*,*) "***** PUT_VAR problem in writediagfi"
366              write(*,*) "***** with ",nom
367              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
368c             call abort
369           endif
370
371          endif !of if (is_master)
372
373!Case of a 2D variable
374!---------------------
375
376        else if (dim.eq.2) then
377
378#ifdef CPP_PARA
379          ! Gather field on a "global" (without redundant longitude) array
380          px2(:)=px(:,1)
381          call Gather(px2,dx2_glop)
382!$OMP MASTER
383          if (is_mpi_master) then
384            call Grid1Dto2D_glo(dx2_glop,dx2_glo)
385            ! copy dx2_glo() to dx2(:) and add redundant longitude
386            dx2(1:iim,:)=dx2_glo(1:iim,:)
387            dx2(iip1,:)=dx2(1,:)
388          endif
389!$OMP END MASTER
390!$OMP BARRIER
391#else
392
393!         Passage variable physique -->  physique dynamique
394!         recast (copy) variable from physics grid to dynamics grid
395
396             DO i=1,iip1
397                dx2(i,1)=px(1,1)
398                dx2(i,jjp1)=px(ngrid,1)
399             ENDDO
400             DO j=2,jjm
401                ig0= 1+(j-2)*iim
402                DO i=1,iim
403                   dx2(i,j)=px(ig0+i,1)
404                ENDDO
405                dx2(iip1,j)=dx2(1,j)
406             ENDDO
407#endif
408
409          if (is_master) then
410           ! only the master writes to output
411!         write (*,*) 'In  writediagfi, on sauve:  ' , nom
412!         write (*,*) 'In  writediagfi. Estimated date = ' ,date
413           ierr= NF_INQ_VARID(nid,nom,varid)
414           if (ierr /= NF_NOERR) then
415! corresponding dimensions
416              ierr= NF_INQ_DIMID(nid,"longitude",id(1))
417              ierr= NF_INQ_DIMID(nid,"latitude",id(2))
418              ierr= NF_INQ_DIMID(nid,"Time",id(3))
419
420! Create the variable if it doesn't exist yet
421
422              write (*,*) "=========================="
423              write (*,*) "DIAGFI: creating variable ",nom
424
425              call def_var(nid,nom,titre,unite,3,id,varid,ierr)
426
427           endif
428
429           corner(1)=1
430           corner(2)=1
431           corner(3)=ntime
432           edges(1)=iip1
433           edges(2)=jjp1
434           edges(3)=1
435
436
437!#ifdef NC_DOUBLE
438!           ierr = NF_PUT_VARA_DOUBLE (nid,varid,corner,edges,dx2)
439!#else         
440           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx2)
441!#endif     
442
443           if (ierr.ne.NF_NOERR) then
444              write(*,*) "***** PUT_VAR matter in writediagfi"
445              write(*,*) "***** with ",nom
446              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
447c             call abort
448           endif
449
450          endif !of if (is_master)
451
452!Case of a 1D variable (ie: a column)
453!---------------------------------------------------
454
455       else if (dim.eq.1) then
456         if (is_parallel) then
457           write(*,*) "writediagfi error: dim=1 not implemented ",
458     &                 "in parallel mode"
459           CALL abort_physiq
460         endif
461!         Passage variable physique -->  physique dynamique
462!         recast (copy) variable from physics grid to dynamics grid
463          do l=1,llm
464            dx1(l)=px(1,l)
465          enddo
466         
467          ierr= NF_INQ_VARID(nid,nom,varid)
468           if (ierr /= NF_NOERR) then
469! corresponding dimensions
470              ierr= NF_INQ_DIMID(nid,"altitude",id(1))
471              ierr= NF_INQ_DIMID(nid,"Time",id(2))
472
473! Create the variable if it doesn't exist yet
474
475              write (*,*) "=========================="
476              write (*,*) "DIAGFI: creating variable ",nom
477
478              call def_var(nid,nom,titre,unite,2,id,varid,ierr)
479             
480           endif
481           
482           corner(1)=1
483           corner(2)=ntime
484           
485           edges(1)=llm
486           edges(2)=1
487!#ifdef NC_DOUBLE
488!           ierr= NF_PUT_VARA_DOUBLE(nid,varid,corner,edges,dx1)
489!#else
490           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx1)
491!#endif
492
493           if (ierr.ne.NF_NOERR) then
494              write(*,*) "***** PUT_VAR problem in writediagfi"
495              write(*,*) "***** with ",nom
496              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
497c             call abort
498           endif 
499
500!Case of a 0D variable (ie: a time-dependent scalar)
501!---------------------------------------------------
502
503        else if (dim.eq.0) then
504
505           dx0 = px (1,1)
506
507          if (is_master) then
508           ! only the master writes to output
509           ierr= NF_INQ_VARID(nid,nom,varid)
510           if (ierr /= NF_NOERR) then
511! corresponding dimensions
512              ierr= NF_INQ_DIMID(nid,"Time",id(1))
513
514! Create the variable if it doesn't exist yet
515
516              write (*,*) "=========================="
517              write (*,*) "DIAGFI: creating variable ",nom
518
519              call def_var(nid,nom,titre,unite,1,id,varid,ierr)
520
521           endif
522
523           corner(1)=ntime
524           edges(1)=1
525
526!#ifdef NC_DOUBLE
527!           ierr = NF_PUT_VARA_DOUBLE (nid,varid,corner,edges,dx0) 
528!#else
529           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx0)
530!#endif
531           if (ierr.ne.NF_NOERR) then
532              write(*,*) "***** PUT_VAR matter in writediagfi"
533              write(*,*) "***** with ",nom
534              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
535c             call abort
536           endif
537
538          endif !of if (is_master)
539
540        endif ! of if (dim.eq.3) elseif(dim.eq.2)...
541
542      endif ! of if ( MOD(zitau+1,irythme) .eq.0.)
543
544      if (is_master) then
545        ierr= NF_CLOSE(nid)
546      endif
547
548#endif
549      end
Note: See TracBrowser for help on using the repository browser.