source: trunk/LMDZ.GENERIC/libf/phystd/writediagfi.F @ 983

Last change on this file since 983 was 969, checked in by emillour, 12 years ago

Generic/Universal? GCM:
Further cleanup in outputs: enable output of a scalar in writediagfi in parallel mode and remove obsolete io* routines.
EM

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