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

Last change on this file since 245 was 245, checked in by ymipsl, 10 years ago
  • One call for initialize physics from dynamico
  • mpi_root renamed into mpi_master due to conflict with an existaing symbol from the mpi library

==> mpi_root => mpi_master, is_mpi_root => is_mpi_master, is_omp_root => is_omp_master

YM

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