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

Last change on this file since 966 was 965, checked in by emillour, 12 years ago

Common dynamics and generic/universal GCM:

  • LMDZ.COMMON: minor bug fix on the computation of physics mesh area in gcm.F
  • LMDZ.UNIVERSAL: missing clean initialization of tab_cntrl(:) array in phyredem.F90
  • LMDZ.GENERIC: minor bug fix in hydrol.F90, only output runoff if it is used. Update output routines so that all outputs files (stats, diagfi.nc, diagsoil.nc, diagspecIR.nc and diagspecVI.nc) can be generated when running LMDZ.UNIVERSAL in MPI mode.

EM

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