source: trunk/LMDZ.MARS/libf/phymars/writediagfi.F @ 320

Last change on this file since 320 was 320, checked in by emillour, 13 years ago

Mars GCM:

  • added FF's upgrade of writediagfi. Now, if at runtime there is a diagfi.def file, it should contain the list of variables (1 per line) than will be put in the diagfi.nc file. If there is no diagfi.def file, then all variables are put in the diagfi.nc file (as was the case before).

EM

File size: 14.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      implicit none
43
44! Commons
45#include "dimensions.h"
46#include "dimphys.h"
47#include "paramet.h"
48#include "control.h"
49#include "comvert.h"
50#include "comgeom.h"
51#include "description.h"
52#include "netcdf.inc"
53#include "temps.h"
54#include "surfdat.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
71      REAL phis(ip1jmp1)
72
73      integer irythme
74      integer ierr,ierr2
75      integer iq
76      integer i,j,l,zmax , ig0
77
78      integer,save :: zitau=0
79      character(len=20),save :: firstnom='1234567890'
80
81! Ajouts
82      integer, save :: ntime=0
83      integer :: idim,varid
84      integer :: nid
85      character(len=*),parameter :: fichnom="diagfi.nc"
86      integer, dimension(4) :: id
87      integer, dimension(4) :: edges,corner
88
89! Added to use diagfi.def to select output variable
90      logical,save :: diagfi_def
91      logical :: getout
92      integer,save :: n_nom_def
93      integer :: n
94      integer,parameter :: n_nom_def_max=99
95      character(len=20),save :: nom_def(n_nom_def_max)
96      logical,save :: firstcall=.true.
97     
98#ifndef MESOSCALE
99!***************************************************************
100!Sortie des variables au rythme voulu
101
102      irythme = int(ecritphy) ! sortie au rythme de ecritphy
103!     irythme = iconser  ! sortie au rythme des variables de controle
104!     irythme = iphysiq  ! sortie a tous les pas physique
105!     irythme = iecri*day_step ! sortie au rythme des fichiers histoires
106!     irythme = periodav*day_step ! sortie au rythme des fichiers histmoy
107
108!***************************************************************
109
110! At very first call, check if there is a "diagfi.def" to use and read it
111! -----------------------------------------------------------------------
112      IF (firstcall) THEN
113         firstcall=.false.
114
115  !      Open diagfi.def definition file if there is one:
116         open(99,file="diagfi.def",status='old',form='formatted',
117     s   iostat=ierr2)
118
119         if(ierr2.eq.0) then
120            diagfi_def=.true.
121            write(*,*) "******************"
122            write(*,*) "Reading diagfi.def"
123            write(*,*) "******************"
124            do n=1,n_nom_def_max
125              read(99,fmt='(a)',end=88) nom_def(n)
126              write(*,*) 'Output in diagfi: ', trim(nom_def(n))
127            end do
128 88         continue
129            if (n.ge.n_nom_def_max) then
130               write(*,*)"n_nom_def_max too small in writediagfi.F:",n
131               stop
132            end if
133            n_nom_def=n-1
134            close(99)
135         else
136            diagfi_def=.false.
137         endif
138      END IF ! of IF (firstcall)
139
140! Get out of write_diagfi if there is diagfi.def AND variable not listed
141!  ---------------------------------------------------------------------
142      if (diagfi_def) then
143          getout=.true.
144          do n=1,n_nom_def
145             if(trim(nom_def(n)).eq.nom) getout=.false.
146          end do
147          if (getout) return
148      end if
149
150! Initialisation of 'firstnom' and create/open the "diagfi.nc" NetCDF file
151! ------------------------------------------------------------------------
152! (at very first call to the subroutine, in accordance with diagfi.def)
153
154      if (firstnom.eq.'1234567890') then ! .true. for the very first valid
155      !   call to this subroutine; now set 'firstnom'
156         firstnom = nom
157         ! just to be sure, check that firstnom is large enough to hold nom
158         if (len_trim(firstnom).lt.len_trim(nom)) then
159           write(*,*) "writediagfi: Error !!!"
160           write(*,*) "   firstnom string not long enough!!"
161           write(*,*) "   increase its size to at least ",len_trim(nom)
162           stop
163         endif
164
165         ! Create the NetCDF file
166         ierr = NF_CREATE(fichnom, NF_CLOBBER, nid)
167         ! Define the 'Time' dimension
168         ierr = nf_def_dim(nid,"Time",NF_UNLIMITED,idim)
169         ! Define the 'Time' variable
170!#ifdef NC_DOUBLE
171!         ierr = NF_DEF_VAR (nid, "Time", NF_DOUBLE, 1, idim,varid)
172!#else
173         ierr = NF_DEF_VAR (nid, "Time", NF_FLOAT, 1, idim,varid)
174!#endif
175         ! Add a long_name attribute
176         ierr = NF_PUT_ATT_TEXT (nid, varid, "long_name",
177     .          4,"Time")
178         ! Add a units attribute
179         ierr = NF_PUT_ATT_TEXT(nid, varid,'units',29,
180     .          "days since 0000-00-0 00:00:00")
181         ! Switch out of NetCDF Define mode
182         ierr = NF_ENDDEF(nid)
183
184         ! write "header" of file (longitudes, latitudes, geopotential, ...)
185         call gr_fi_dyn(1,ngrid,iip1,jjp1,phisfi,phis)
186         call iniwrite(nid,day_ini,phis)
187         
188         zitau = -1 ! initialize zitau
189      else
190         ! Open the NetCDF file
191         ierr = NF_OPEN(fichnom,NF_WRITE,nid)
192      endif ! if (firstnom.eq.'1234567890')
193
194      if (ngridmx.eq.1) then
195        ! in testphys1d, for the 1d version of the GCM, iphysiq and irythme
196        ! are undefined; so set them to 1
197        iphysiq=1
198        irythme=1
199        ! NB:
200      endif
201
202! Increment time index 'zitau' if it is the "fist call" (at given time level)
203! to writediagfi
204!------------------------------------------------------------------------
205      if (nom.eq.firstnom) then
206          zitau = zitau + iphysiq
207      end if
208
209!--------------------------------------------------------
210! Write the variables to output file if it's time to do so
211!--------------------------------------------------------
212
213      if ( MOD(zitau+1,irythme) .eq.0.) then
214
215! Compute/write/extend 'Time' coordinate (date given in days)
216! (done every "first call" (at given time level) to writediagfi)
217! Note: date is incremented as 1 step ahead of physics time
218!--------------------------------------------------------
219
220        if (nom.eq.firstnom) then
221        ! We have identified a "first call" (at given date)
222           ntime=ntime+1 ! increment # of stored time steps
223           ! compute corresponding date (in days and fractions thereof)
224           date= float (zitau +1)/float (day_step)
225           ! Get NetCDF ID of 'Time' variable
226           ierr= NF_INQ_VARID(nid,"Time",varid)
227           ! Write (append) the new date to the 'Time' array
228!#ifdef NC_DOUBLE
229!           ierr= NF_PUT_VARA_DOUBLE(nid,varid,ntime,1,date)
230!#else
231           ierr= NF_PUT_VARA_REAL(nid,varid,ntime,1,date)
232!#endif
233           if (ierr.ne.NF_NOERR) then
234              write(*,*) "***** PUT_VAR matter in writediagfi_nc"
235              write(*,*) "***** with time"
236              write(*,*) 'ierr=', ierr   
237c             call abort
238           endif
239
240           write(6,*)'WRITEDIAGFI: date= ', date
241        end if ! of if (nom.eq.firstnom)
242
243!Case of a 3D variable
244!---------------------
245        if (dim.eq.3) then
246
247!         Passage variable physique -->  variable dynamique
248!         recast (copy) variable from physics grid to dynamics grid
249           DO l=1,llm
250             DO i=1,iip1
251                dx3(i,1,l)=px(1,l)
252                dx3(i,jjp1,l)=px(ngrid,l)
253             ENDDO
254             DO j=2,jjm
255                ig0= 1+(j-2)*iim
256                DO i=1,iim
257                   dx3(i,j,l)=px(ig0+i,l)
258                ENDDO
259                dx3(iip1,j,l)=dx3(1,j,l)
260             ENDDO
261           ENDDO
262
263!         Ecriture du champs
264
265!         write (*,*) 'In  writediagfi, on sauve:  ' , nom
266!         write (*,*) 'In  writediagfi. Estimated date = ' ,date
267! name of the variable
268           ierr= NF_INQ_VARID(nid,nom,varid)
269           if (ierr /= NF_NOERR) then
270! corresponding dimensions
271              ierr= NF_INQ_DIMID(nid,"longitude",id(1))
272              ierr= NF_INQ_DIMID(nid,"latitude",id(2))
273              ierr= NF_INQ_DIMID(nid,"altitude",id(3))
274              ierr= NF_INQ_DIMID(nid,"Time",id(4))
275
276! Create the variable if it doesn't exist yet
277
278              write (*,*) "=========================="
279              write (*,*) "DIAGFI: creating variable ",nom
280              call def_var(nid,nom,titre,unite,4,id,varid,ierr)
281
282           endif
283
284           corner(1)=1
285           corner(2)=1
286           corner(3)=1
287           corner(4)=ntime
288
289           edges(1)=iip1
290           edges(2)=jjp1
291           edges(3)=llm
292           edges(4)=1
293!#ifdef NC_DOUBLE
294!           ierr= NF_PUT_VARA_DOUBLE(nid,varid,corner,edges,dx3)
295!#else
296           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx3)
297!#endif
298
299           if (ierr.ne.NF_NOERR) then
300              write(*,*) "***** PUT_VAR problem in writediagfi"
301              write(*,*) "***** with ",nom
302              write(*,*) 'ierr=', ierr
303c             call abort
304           endif
305
306!Case of a 2D variable
307!---------------------
308
309        else if (dim.eq.2) then
310
311!         Passage variable physique -->  physique dynamique
312!         recast (copy) variable from physics grid to dynamics grid
313
314             DO i=1,iip1
315                dx2(i,1)=px(1,1)
316                dx2(i,jjp1)=px(ngrid,1)
317             ENDDO
318             DO j=2,jjm
319                ig0= 1+(j-2)*iim
320                DO i=1,iim
321                   dx2(i,j)=px(ig0+i,1)
322                ENDDO
323                dx2(iip1,j)=dx2(1,j)
324             ENDDO
325
326!         write (*,*) 'In  writediagfi, on sauve:  ' , nom
327!         write (*,*) 'In  writediagfi. Estimated date = ' ,date
328           ierr= NF_INQ_VARID(nid,nom,varid)
329           if (ierr /= NF_NOERR) then
330! corresponding dimensions
331              ierr= NF_INQ_DIMID(nid,"longitude",id(1))
332              ierr= NF_INQ_DIMID(nid,"latitude",id(2))
333              ierr= NF_INQ_DIMID(nid,"Time",id(3))
334
335! Create the variable if it doesn't exist yet
336
337              write (*,*) "=========================="
338              write (*,*) "DIAGFI: creating variable ",nom
339
340              call def_var(nid,nom,titre,unite,3,id,varid,ierr)
341
342           endif
343
344           corner(1)=1
345           corner(2)=1
346           corner(3)=ntime
347           edges(1)=iip1
348           edges(2)=jjp1
349           edges(3)=1
350
351
352!#ifdef NC_DOUBLE
353!           ierr = NF_PUT_VARA_DOUBLE (nid,varid,corner,edges,dx2)
354!#else         
355           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx2)
356!#endif     
357
358           if (ierr.ne.NF_NOERR) then
359              write(*,*) "***** PUT_VAR matter in writediagfi"
360              write(*,*) "***** with ",nom
361              write(*,*) 'ierr=', ierr
362c             call abort
363           endif
364
365!Case of a 1D variable (ie: a column)
366!---------------------------------------------------
367
368       else if (dim.eq.1) then
369!         Passage variable physique -->  physique dynamique
370!         recast (copy) variable from physics grid to dynamics grid
371          do l=1,llm
372            dx1(l)=px(1,l)
373          enddo
374         
375          ierr= NF_INQ_VARID(nid,nom,varid)
376           if (ierr /= NF_NOERR) then
377! corresponding dimensions
378              ierr= NF_INQ_DIMID(nid,"altitude",id(1))
379              ierr= NF_INQ_DIMID(nid,"Time",id(2))
380
381! Create the variable if it doesn't exist yet
382
383              write (*,*) "=========================="
384              write (*,*) "DIAGFI: creating variable ",nom
385
386              call def_var(nid,nom,titre,unite,2,id,varid,ierr)
387             
388           endif
389           
390           corner(1)=1
391           corner(2)=ntime
392           
393           edges(1)=llm
394           edges(2)=1
395!#ifdef NC_DOUBLE
396!           ierr= NF_PUT_VARA_DOUBLE(nid,varid,corner,edges,dx1)
397!#else
398           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx1)
399!#endif
400
401           if (ierr.ne.NF_NOERR) then
402              write(*,*) "***** PUT_VAR problem in writediagfi"
403              write(*,*) "***** with ",nom
404              write(*,*) 'ierr=', ierr
405c             call abort
406           endif
407
408!Case of a 0D variable (ie: a time-dependent scalar)
409!---------------------------------------------------
410
411        else if (dim.eq.0) then
412           dx0 = px (1,1)
413
414           ierr= NF_INQ_VARID(nid,nom,varid)
415           if (ierr /= NF_NOERR) then
416! corresponding dimensions
417              ierr= NF_INQ_DIMID(nid,"Time",id(1))
418
419! Create the variable if it doesn't exist yet
420
421              write (*,*) "=========================="
422              write (*,*) "DIAGFI: creating variable ",nom
423
424              call def_var(nid,nom,titre,unite,1,id,varid,ierr)
425
426           endif
427
428           corner(1)=ntime
429           edges(1)=1
430
431!#ifdef NC_DOUBLE
432!           ierr = NF_PUT_VARA_DOUBLE (nid,varid,corner,edges,dx0) 
433!#else
434           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx0)
435!#endif
436           if (ierr.ne.NF_NOERR) then
437              write(*,*) "***** PUT_VAR matter in writediagfi"
438              write(*,*) "***** with ",nom
439              write(*,*) 'ierr=', ierr
440c             call abort
441           endif
442
443        endif ! of if (dim.eq.3) elseif(dim.eq.2)...
444
445      endif ! of if ( MOD(zitau+1,irythme) .eq.0.)
446
447      ierr= NF_CLOSE(nid)
448
449#endif
450      end
Note: See TracBrowser for help on using the repository browser.