source: trunk/LMDZ.PLUTO.old/libf/phypluto/writediagfi.F @ 3436

Last change on this file since 3436 was 3175, checked in by emillour, 11 months ago

Pluto PCM:
Add the old Pluto LMDZ for reference (required prior step to making
an LMDZ.PLUTO using the same framework as the other physics packages).
TB+EM

File size: 15.0 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      REAL phis(iip1,jjp1)   !TB18c
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
82! Ajouts
83      integer, save :: ntime=0
84      integer :: idim,varid
85      integer :: nid
86      character(len=*),parameter :: fichnom="diagfi.nc"
87      integer, dimension(4) :: id
88      integer, dimension(4) :: edges,corner
89
90! Added to use diagfi.def to select output variable
91      logical,save :: diagfi_def
92      logical :: getout
93      integer,save :: n_nom_def
94      integer :: n
95      integer,parameter :: n_nom_def_max=99
96      character(len=20),save :: nom_def(n_nom_def_max)
97      logical,save :: firstcall=.true.
98     
99#ifndef MESOSCALE
100!***************************************************************
101!Sortie des variables au rythme voulu
102
103      irythme = int(ecritphy) ! sortie au rythme de ecritphy
104!     irythme = iconser  ! sortie au rythme des variables de controle
105!     irythme = iphysiq  ! sortie a tous les pas physique
106!     irythme = iecri*day_step ! sortie au rythme des fichiers histoires
107!     irythme = periodav*day_step ! sortie au rythme des fichiers histmoy
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         ! TB22 write "header" of file (longitudes, latitudes, geopotential, ...)
185         !if (ngrid.eq.1) then
186         call gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
187         !else
188         !  call gr_fi_dyn2d(1,ngridmx,iip1,jjp1,phisfi,phis)
189         !endif
190         call iniwrite(nid,day_ini,phis)
191         
192         zitau = -1 ! initialize zitau
193      else
194         ! Open the NetCDF file
195         ierr = NF_OPEN(fichnom,NF_WRITE,nid)
196      endif ! if (firstnom.eq.'1234567890')
197
198      if (ngridmx.eq.1) then
199        ! in testphys1d, for the 1d version of the GCM, iphysiq and irythme
200        ! are undefined; so set them to 1
201        iphysiq=1
202        irythme=1
203        ! NB:
204      endif
205
206! Increment time index 'zitau' if it is the "fist call" (at given time level)
207! to writediagfi
208!------------------------------------------------------------------------
209      if (nom.eq.firstnom) then
210          zitau = zitau + iphysiq
211      end if
212
213!--------------------------------------------------------
214! Write the variables to output file if it's time to do so
215!--------------------------------------------------------
216
217      if ( MOD(zitau+1,irythme) .eq.0.) then
218
219! Compute/write/extend 'Time' coordinate (date given in days)
220! (done every "first call" (at given time level) to writediagfi)
221! Note: date is incremented as 1 step ahead of physics time
222!--------------------------------------------------------
223
224        if (nom.eq.firstnom) then
225        ! We have identified a "first call" (at given date)
226           ntime=ntime+1 ! increment # of stored time steps
227           ! compute corresponding date (in days and fractions thereof)
228           date= float (zitau +1)/float (day_step)
229           ! Get NetCDF ID of 'Time' variable
230           ierr= NF_INQ_VARID(nid,"Time",varid)
231           ! Write (append) the new date to the 'Time' array
232!#ifdef NC_DOUBLE
233!           ierr= NF_PUT_VARA_DOUBLE(nid,varid,ntime,1,date)
234!#else
235           ierr= NF_PUT_VARA_REAL(nid,varid,ntime,1,date)
236!#endif
237           if (ierr.ne.NF_NOERR) then
238              write(*,*) "***** PUT_VAR matter in writediagfi_nc"
239              write(*,*) "***** with time"
240              write(*,*) 'ierr4=', ierr,nid,varid,ntime,date   
241              call abort
242           endif
243
244           write(6,*)'WRITEDIAGFI: date= ', date
245        end if ! of if (nom.eq.firstnom)
246
247!Case of a 3D variable
248!---------------------
249        if (dim.eq.3) then
250
251!         Passage variable physique -->  variable dynamique
252!         recast (copy) variable from physics grid to dynamics grid
253           DO l=1,llm
254             DO i=1,iip1
255                dx3(i,1,l)=px(1,l)
256                dx3(i,jjp1,l)=px(ngrid,l)
257             ENDDO
258             DO j=2,jjm
259                ig0= 1+(j-2)*iim
260                DO i=1,iim
261                   dx3(i,j,l)=px(ig0+i,l)
262                ENDDO
263                dx3(iip1,j,l)=dx3(1,j,l)
264             ENDDO
265           ENDDO
266
267!         Ecriture du champs
268
269!         write (*,*) 'In  writediagfi, on sauve:  ' , nom
270!         write (*,*) 'In  writediagfi. Estimated date = ' ,date
271! name of the variable
272           ierr= NF_INQ_VARID(nid,nom,varid)
273           if (ierr /= NF_NOERR) then
274! corresponding dimensions
275              ierr= NF_INQ_DIMID(nid,"longitude",id(1))
276              ierr= NF_INQ_DIMID(nid,"latitude",id(2))
277              ierr= NF_INQ_DIMID(nid,"altitude",id(3))
278              ierr= NF_INQ_DIMID(nid,"Time",id(4))
279
280! Create the variable if it doesn't exist yet
281
282              write (*,*) "=========================="
283              write (*,*) "DIAGFI: creating variable ",nom
284              call def_var(nid,nom,titre,unite,4,id,varid,ierr)
285
286           endif
287
288           corner(1)=1
289           corner(2)=1
290           corner(3)=1
291           corner(4)=ntime
292
293           edges(1)=iip1
294           edges(2)=jjp1
295           edges(3)=llm
296           edges(4)=1
297!#ifdef NC_DOUBLE
298!           ierr= NF_PUT_VARA_DOUBLE(nid,varid,corner,edges,dx3)
299!#else
300           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx3)
301!#endif
302
303           if (ierr.ne.NF_NOERR) then
304              write(*,*) "***** PUT_VAR problem in writediagfi"
305              write(*,*) "***** with ",nom
306              write(*,*) 'ierr5=', ierr
307              call abort
308           endif
309
310!Case of a 2D variable
311!---------------------
312
313        else if (dim.eq.2) then
314
315!         Passage variable physique -->  physique dynamique
316!         recast (copy) variable from physics grid to dynamics grid
317
318             DO i=1,iip1
319                dx2(i,1)=px(1,1)
320                dx2(i,jjp1)=px(ngrid,1)
321             ENDDO
322             DO j=2,jjm
323                ig0= 1+(j-2)*iim
324                DO i=1,iim
325                   dx2(i,j)=px(ig0+i,1)
326                ENDDO
327                dx2(iip1,j)=dx2(1,j)
328             ENDDO
329
330!         write (*,*) 'In  writediagfi, on sauve:  ' , nom
331!         write (*,*) 'In  writediagfi. Estimated date = ' ,date
332           ierr= NF_INQ_VARID(nid,nom,varid)
333           if (ierr /= NF_NOERR) then
334! corresponding dimensions
335              ierr= NF_INQ_DIMID(nid,"longitude",id(1))
336              ierr= NF_INQ_DIMID(nid,"latitude",id(2))
337              ierr= NF_INQ_DIMID(nid,"Time",id(3))
338
339! Create the variable if it doesn't exist yet
340
341              write (*,*) "=========================="
342              write (*,*) "DIAGFI: creating variable ",nom
343
344              call def_var(nid,nom,titre,unite,3,id,varid,ierr)
345
346           endif
347
348           corner(1)=1
349           corner(2)=1
350           corner(3)=ntime
351           edges(1)=iip1
352           edges(2)=jjp1
353           edges(3)=1
354
355
356!#ifdef NC_DOUBLE
357!           ierr = NF_PUT_VARA_DOUBLE (nid,varid,corner,edges,dx2)
358!#else         
359           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx2)
360!#endif     
361
362           if (ierr.ne.NF_NOERR) then
363              write(*,*) "***** PUT_VAR matter in writediagfi"
364              write(*,*) "***** with ",nom
365              write(*,*) 'ierr1=', ierr,nid,varid,corner,edges,dx2
366              call abort
367           endif
368
369!Case of a 1D variable (ie: a column)
370!---------------------------------------------------
371
372       else if (dim.eq.1) then
373!         Passage variable physique -->  physique dynamique
374!         recast (copy) variable from physics grid to dynamics grid
375          do l=1,llm
376            dx1(l)=px(1,l)
377          enddo
378         
379          ierr= NF_INQ_VARID(nid,nom,varid)
380           if (ierr /= NF_NOERR) then
381! corresponding dimensions
382              ierr= NF_INQ_DIMID(nid,"altitude",id(1))
383              ierr= NF_INQ_DIMID(nid,"Time",id(2))
384
385! Create the variable if it doesn't exist yet
386
387              write (*,*) "=========================="
388              write (*,*) "DIAGFI: creating variable ",nom
389
390              call def_var(nid,nom,titre,unite,2,id,varid,ierr)
391             
392           endif
393           
394           corner(1)=1
395           corner(2)=ntime
396           
397           edges(1)=llm
398           edges(2)=1
399!#ifdef NC_DOUBLE
400!           ierr= NF_PUT_VARA_DOUBLE(nid,varid,corner,edges,dx1)
401!#else
402           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx1)
403!#endif
404
405           if (ierr.ne.NF_NOERR) then
406              write(*,*) "***** PUT_VAR problem in writediagfi"
407              write(*,*) "***** with ",nom
408              write(*,*) 'ierr2=', ierr,nid,corner,edges,dx1
409              call abort
410           endif
411
412!Case of a 0D variable (ie: a time-dependent scalar)
413!---------------------------------------------------
414
415        else if (dim.eq.0) then
416           !write (*,*) "TB22 w2 ",nom,px(1,1)
417           !write (*,*) "TB22 px ",px(1,1)
418           dx0 = px(1,1)
419
420           ierr= NF_INQ_VARID(nid,nom,varid)
421           
422           if (ierr /= NF_NOERR) then
423! corresponding dimensions
424              ierr= NF_INQ_DIMID(nid,"Time",id(1))
425
426! Create the variable if it doesn't exist yet
427
428              write (*,*) "=========================="
429              write (*,*) "DIAGFI: creating variable ",nom
430
431              call def_var(nid,nom,titre,unite,1,id,varid,ierr)
432
433           endif
434
435           corner(1)=ntime
436           edges(1)=1
437
438!#ifdef NC_DOUBLE
439!           ierr = NF_PUT_VARA_DOUBLE (nid,varid,corner,edges,dx0) 
440!#else
441           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx0)
442!#endif
443           if (ierr.ne.NF_NOERR) then
444              write(*,*) "***** PUT_VAR matter in writediagfi"
445              write(*,*) "***** with ",nom
446              write(*,*) 'ierr3=', ierr,nid,varid,corner,edges,dx0
447              call abort
448           endif
449
450        endif ! of if (dim.eq.3) elseif(dim.eq.2)...
451
452      endif ! of if ( MOD(zitau+1,irythme) .eq.0.)
453
454      ierr= NF_CLOSE(nid)
455
456#endif
457      end
Note: See TracBrowser for help on using the repository browser.