source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/config/ppsrc/phys/writediagfi.f @ 224

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