source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/writediagfi.F @ 224

Last change on this file since 224 was 222, checked in by ymipsl, 10 years ago

Creating temporary dynamico/lmdz/saturn branche

YM

File size: 17.5 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#ifndef MESOSCALE
99
100#ifdef CPP_PARA
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#else
111      real phisfi_glo(ngrid) ! surface geopotential on global physics grid
112#endif
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#ifdef CPP_PARA
183          ! Gather phisfi() geopotential on physics grid
184          call Gather(phisfi,phisfi_glo)
185#else
186         phisfi_glo(:)=phisfi(:)
187#endif
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 NC_DOUBLE
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 NC_DOUBLE
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#ifdef CPP_PARA
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#else
300!         Passage variable physique -->  variable dynamique
301!         recast (copy) variable from physics grid to dynamics grid
302           DO l=1,llm
303             DO i=1,iip1
304                dx3(i,1,l)=px(1,l)
305                dx3(i,jjp1,l)=px(ngrid,l)
306             ENDDO
307             DO j=2,jjm
308                ig0= 1+(j-2)*iim
309                DO i=1,iim
310                   dx3(i,j,l)=px(ig0+i,l)
311                ENDDO
312                dx3(iip1,j,l)=dx3(1,j,l)
313             ENDDO
314           ENDDO
315#endif
316!         Ecriture du champs
317
318          if (is_master) then
319           ! only the master writes to output
320! name of the variable
321           ierr= NF_INQ_VARID(nid,nom,varid)
322           if (ierr /= NF_NOERR) then
323! corresponding dimensions
324              ierr= NF_INQ_DIMID(nid,"longitude",id(1))
325              ierr= NF_INQ_DIMID(nid,"latitude",id(2))
326              ierr= NF_INQ_DIMID(nid,"altitude",id(3))
327              ierr= NF_INQ_DIMID(nid,"Time",id(4))
328
329! Create the variable if it doesn't exist yet
330
331              write (*,*) "=========================="
332              write (*,*) "DIAGFI: creating variable ",nom
333              call def_var(nid,nom,titre,unite,4,id,varid,ierr)
334
335           endif
336
337           corner(1)=1
338           corner(2)=1
339           corner(3)=1
340           corner(4)=ntime
341
342           edges(1)=iip1
343           edges(2)=jjp1
344           edges(3)=llm
345           edges(4)=1
346!#ifdef NC_DOUBLE
347!           ierr= NF_PUT_VARA_DOUBLE(nid,varid,corner,edges,dx3)
348!#else
349!           write(*,*)"test:  nid=",nid," varid=",varid
350!           write(*,*)"       corner()=",corner
351!           write(*,*)"       edges()=",edges
352!           write(*,*)"       dx3()=",dx3
353           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx3)
354!#endif
355
356           if (ierr.ne.NF_NOERR) then
357              write(*,*) "***** PUT_VAR problem in writediagfi"
358              write(*,*) "***** with ",nom
359              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
360c             call abort
361           endif
362
363          endif !of if (is_master)
364
365!Case of a 2D variable
366!---------------------
367
368        else if (dim.eq.2) then
369
370#ifdef CPP_PARA
371          ! Gather field on a "global" (without redundant longitude) array
372          px2(:)=px(:,1)
373          call Gather(px2,dx2_glop)
374!$OMP MASTER
375          if (is_mpi_root) then
376            call Grid1Dto2D_glo(dx2_glop,dx2_glo)
377            ! copy dx2_glo() to dx2(:) and add redundant longitude
378            dx2(1:iim,:)=dx2_glo(1:iim,:)
379            dx2(iip1,:)=dx2(1,:)
380          endif
381!$OMP END MASTER
382!$OMP BARRIER
383#else
384
385!         Passage variable physique -->  physique dynamique
386!         recast (copy) variable from physics grid to dynamics grid
387
388             DO i=1,iip1
389                dx2(i,1)=px(1,1)
390                dx2(i,jjp1)=px(ngrid,1)
391             ENDDO
392             DO j=2,jjm
393                ig0= 1+(j-2)*iim
394                DO i=1,iim
395                   dx2(i,j)=px(ig0+i,1)
396                ENDDO
397                dx2(iip1,j)=dx2(1,j)
398             ENDDO
399#endif
400
401          if (is_master) then
402           ! only the master writes to output
403!         write (*,*) 'In  writediagfi, on sauve:  ' , nom
404!         write (*,*) 'In  writediagfi. Estimated date = ' ,date
405           ierr= NF_INQ_VARID(nid,nom,varid)
406           if (ierr /= NF_NOERR) then
407! corresponding dimensions
408              ierr= NF_INQ_DIMID(nid,"longitude",id(1))
409              ierr= NF_INQ_DIMID(nid,"latitude",id(2))
410              ierr= NF_INQ_DIMID(nid,"Time",id(3))
411
412! Create the variable if it doesn't exist yet
413
414              write (*,*) "=========================="
415              write (*,*) "DIAGFI: creating variable ",nom
416
417              call def_var(nid,nom,titre,unite,3,id,varid,ierr)
418
419           endif
420
421           corner(1)=1
422           corner(2)=1
423           corner(3)=ntime
424           edges(1)=iip1
425           edges(2)=jjp1
426           edges(3)=1
427
428
429!#ifdef NC_DOUBLE
430!           ierr = NF_PUT_VARA_DOUBLE (nid,varid,corner,edges,dx2)
431!#else         
432           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx2)
433!#endif     
434
435           if (ierr.ne.NF_NOERR) then
436              write(*,*) "***** PUT_VAR matter in writediagfi"
437              write(*,*) "***** with ",nom
438              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
439c             call abort
440           endif
441
442          endif !of if (is_master)
443
444!Case of a 1D variable (ie: a column)
445!---------------------------------------------------
446
447       else if (dim.eq.1) then
448         if (is_parallel) then
449           write(*,*) "writediagfi error: dim=1 not implemented ",
450     &                 "in parallel mode"
451           stop
452         endif
453!         Passage variable physique -->  physique dynamique
454!         recast (copy) variable from physics grid to dynamics grid
455          do l=1,llm
456            dx1(l)=px(1,l)
457          enddo
458         
459          ierr= NF_INQ_VARID(nid,nom,varid)
460           if (ierr /= NF_NOERR) then
461! corresponding dimensions
462              ierr= NF_INQ_DIMID(nid,"altitude",id(1))
463              ierr= NF_INQ_DIMID(nid,"Time",id(2))
464
465! Create the variable if it doesn't exist yet
466
467              write (*,*) "=========================="
468              write (*,*) "DIAGFI: creating variable ",nom
469
470              call def_var(nid,nom,titre,unite,2,id,varid,ierr)
471             
472           endif
473           
474           corner(1)=1
475           corner(2)=ntime
476           
477           edges(1)=llm
478           edges(2)=1
479!#ifdef NC_DOUBLE
480!           ierr= NF_PUT_VARA_DOUBLE(nid,varid,corner,edges,dx1)
481!#else
482           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx1)
483!#endif
484
485           if (ierr.ne.NF_NOERR) then
486              write(*,*) "***** PUT_VAR problem in writediagfi"
487              write(*,*) "***** with ",nom
488              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
489c             call abort
490           endif 
491
492!Case of a 0D variable (ie: a time-dependent scalar)
493!---------------------------------------------------
494
495        else if (dim.eq.0) then
496
497           dx0 = px (1,1)
498
499          if (is_master) then
500           ! only the master writes to output
501           ierr= NF_INQ_VARID(nid,nom,varid)
502           if (ierr /= NF_NOERR) then
503! corresponding dimensions
504              ierr= NF_INQ_DIMID(nid,"Time",id(1))
505
506! Create the variable if it doesn't exist yet
507
508              write (*,*) "=========================="
509              write (*,*) "DIAGFI: creating variable ",nom
510
511              call def_var(nid,nom,titre,unite,1,id,varid,ierr)
512
513           endif
514
515           corner(1)=ntime
516           edges(1)=1
517
518!#ifdef NC_DOUBLE
519!           ierr = NF_PUT_VARA_DOUBLE (nid,varid,corner,edges,dx0) 
520!#else
521           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx0)
522!#endif
523           if (ierr.ne.NF_NOERR) then
524              write(*,*) "***** PUT_VAR matter in writediagfi"
525              write(*,*) "***** with ",nom
526              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
527c             call abort
528           endif
529
530          endif !of if (is_master)
531
532        endif ! of if (dim.eq.3) elseif(dim.eq.2)...
533
534      endif ! of if ( MOD(zitau+1,irythme) .eq.0.)
535
536      if (is_master) then
537        ierr= NF_CLOSE(nid)
538      endif
539
540#endif
541      end
Note: See TracBrowser for help on using the repository browser.