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

Last change on this file since 228 was 228, checked in by aslmd, 13 years ago

MESOSCALE: corrected the diagfi bug. updated test case and check those are still working OK. other changes are cosmetic

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