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

Last change on this file since 988 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
RevLine 
[135]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)
[323]8!  (ou encore 1d, dans le cas de testphys1d, pour sortir une colonne)
[135]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
[323]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
[135]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!=================================================================
[787]41
42      USE surfdat_h
[965]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
[135]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:
[323]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)
[135]66
67! Local variables:
68
[323]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
[135]73
[323]74      real*4,save :: date
[135]75
[965]76      REAL phis(ip1jmp1) ! surface geopotential on extended lonxlat grid
[135]77
78      integer irythme
[323]79      integer ierr,ierr2
[135]80      integer iq
81      integer i,j,l,zmax , ig0
82
[323]83      integer,save :: zitau=0
84      character(len=20),save :: firstnom='1234567890'
[135]85
86! Ajouts
87      integer, save :: ntime=0
88      integer :: idim,varid
89      integer :: nid
[323]90      character(len=*),parameter :: fichnom="diagfi.nc"
[135]91      integer, dimension(4) :: id
92      integer, dimension(4) :: edges,corner
[323]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.
[787]102
[862]103      integer dimvert
104
[965]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
[135]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
[787]132
[323]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.
[135]137
[787]138!      Open diagfi.def definition file if there is one:
[323]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
[862]173
[135]174! Initialisation of 'firstnom' and create/open the "diagfi.nc" NetCDF file
175! ------------------------------------------------------------------------
[323]176! (at very first call to the subroutine, in accordance with diagfi.def)
[135]177
[323]178      if (firstnom.eq.'1234567890') then ! .true. for the very first valid
179      !   call to this subroutine; now set 'firstnom'
[135]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
[862]189         zitau = -1 ! initialize zitau
190
191#ifdef CPP_PARA
[965]192          ! Gather phisfi() geopotential on physics grid
193          call Gather(phisfi,phisfi_glo)
194#else
195         phisfi_glo(:)=phisfi(:)
196#endif
197
[862]198         !! parallel: we cannot use the usual writediagfi method
[965]199!!         call iophys_ini
200         if (is_master) then
201         ! only the master is required to do this
202
[135]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
[323]208!#ifdef NC_DOUBLE
209!         ierr = NF_DEF_VAR (nid, "Time", NF_DOUBLE, 1, idim,varid)
210!#else
[135]211         ierr = NF_DEF_VAR (nid, "Time", NF_FLOAT, 1, idim,varid)
[323]212!#endif
[135]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, ...)
[965]223         call gr_fi_dyn(1,size(phisfi_glo),iip1,jjp1,phisfi_glo,phis)
[135]224         call iniwrite(nid,day_ini,phis)
[862]225
[965]226         endif ! of if (is_master)
227
[135]228      else
[965]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
[135]237      endif ! if (firstnom.eq.'1234567890')
238
[787]239      if (ngrid.eq.1) then
[323]240        ! in testphys1d, for the 1d version of the GCM, iphysiq and irythme
[135]241        ! are undefined; so set them to 1
242        iphysiq=1
[526]243        !JL11 irythme=1
[135]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
[965]265        if (is_master) then
266           ! only the master is required to do this
[135]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
[323]275!#ifdef NC_DOUBLE
276!           ierr= NF_PUT_VARA_DOUBLE(nid,varid,ntime,1,date)
277!#else
[135]278           ierr= NF_PUT_VARA_REAL(nid,varid,ntime,1,date)
[323]279!#endif
[135]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
[965]290        endif ! of if (is_master)
291
[135]292!Case of a 3D variable
293!---------------------
294        if (dim.eq.3) then
295
[965]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
[135]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
[965]324#endif
[135]325!         Ecriture du champs
326
[965]327          if (is_master) then
328           ! only the master writes to output
[135]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
[323]355!#ifdef NC_DOUBLE
356!           ierr= NF_PUT_VARA_DOUBLE(nid,varid,corner,edges,dx3)
357!#else
[965]358!           write(*,*)"test:  nid=",nid," varid=",varid
359!           write(*,*)"       corner()=",corner
360!           write(*,*)"       edges()=",edges
361!           write(*,*)"       dx3()=",dx3
[135]362           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx3)
[323]363!#endif
[135]364
365           if (ierr.ne.NF_NOERR) then
366              write(*,*) "***** PUT_VAR problem in writediagfi"
367              write(*,*) "***** with ",nom
[965]368              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
[135]369c             call abort
370           endif
371
[965]372          endif !of if (is_master)
373
[135]374!Case of a 2D variable
375!---------------------
376
377        else if (dim.eq.2) then
378
[965]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
[135]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
[965]408#endif
[135]409
[965]410          if (is_master) then
411           ! only the master writes to output
[135]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
[323]438!#ifdef NC_DOUBLE
439!           ierr = NF_PUT_VARA_DOUBLE (nid,varid,corner,edges,dx2)
440!#else         
[135]441           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx2)
[323]442!#endif     
[135]443
444           if (ierr.ne.NF_NOERR) then
445              write(*,*) "***** PUT_VAR matter in writediagfi"
446              write(*,*) "***** with ",nom
[965]447              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
[135]448c             call abort
449           endif
450
[965]451          endif !of if (is_master)
452
[135]453!Case of a 1D variable (ie: a column)
454!---------------------------------------------------
455
456       else if (dim.eq.1) then
[965]457         if (is_parallel) then
458           write(*,*) "writediagfi error: dim=1 not implemented ",
459     &                 "in parallel mode"
460           stop
461         endif
[135]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
[323]488!#ifdef NC_DOUBLE
489!           ierr= NF_PUT_VARA_DOUBLE(nid,varid,corner,edges,dx1)
490!#else
[135]491           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx1)
[323]492!#endif
[135]493
494           if (ierr.ne.NF_NOERR) then
495              write(*,*) "***** PUT_VAR problem in writediagfi"
496              write(*,*) "***** with ",nom
[969]497              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
[135]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
[965]505
[135]506           dx0 = px (1,1)
507
[969]508          if (is_master) then
509           ! only the master writes to output
[135]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
[323]527!#ifdef NC_DOUBLE
528!           ierr = NF_PUT_VARA_DOUBLE (nid,varid,corner,edges,dx0) 
529!#else
[135]530           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx0)
[323]531!#endif
[135]532           if (ierr.ne.NF_NOERR) then
533              write(*,*) "***** PUT_VAR matter in writediagfi"
534              write(*,*) "***** with ",nom
[969]535              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
[135]536c             call abort
537           endif
538
[969]539          endif !of if (is_master)
540
[135]541        endif ! of if (dim.eq.3) elseif(dim.eq.2)...
542
543      endif ! of if ( MOD(zitau+1,irythme) .eq.0.)
544
[965]545      if (is_master) then
546        ierr= NF_CLOSE(nid)
547      endif
[787]548
[135]549      end
[965]550
Note: See TracBrowser for help on using the repository browser.