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

Last change on this file since 937 was 862, checked in by aslmd, 12 years ago

LMDZ.UNIVERSAL
LMDZ.GENERIC

Added calls to Frederic Hourdin's subroutines to output physical fields in parallel. This is under precompiling flags CPP_PARA.
This allows for LMDZ.UNIVERSAL users to use writediagfi with parallel computations.
These lines are not compiled by casual users of LMDZ.GENERIC (or users of LMDZ.UNIVERSAL in sequential mode).

The precompiling flag NOWRITEDIAGFI is obsolete and has been deleted.

NOTES:

  • A better cleaner method might be proposed later
  • Added subroutines

A 0 LMDZ.UNIVERSAL/libf/phygeneric/iotd_ecrit.F90
A 0 LMDZ.UNIVERSAL/libf/phygeneric/iophys.F90
A 0 LMDZ.UNIVERSAL/libf/phygeneric/iotd.h
A 0 LMDZ.UNIVERSAL/libf/phygeneric/iotd_ini.F90
A 0 LMDZ.UNIVERSAL/libf/phygeneric/iotd_fin.F90
The iotd_* subroutines are actually supposed to be put in bibio in LMDZ.COMMON
This will be done later once agreed in the team

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