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

Last change on this file since 1477 was 1422, checked in by milmd, 10 years ago

In GENERIC, MARS and COMMON models replace some include files by modules (usefull for decoupling physics with dynamics).

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