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

Last change on this file since 1828 was 1828, checked in by mlefevre, 7 years ago

corrected a weird back-to-the-future bug in the previous commit

File size: 19.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 geometry_mod, only: cell_area
43      use time_phylmdz_mod, only: ecritphy, day_step, iphysiq, day_ini
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     &                              nbp_lon, nbp_lat, nbp_lev
48      implicit none
49
50! Commons
51      include "netcdf.inc"
52
53! Arguments on input:
54      integer,intent(in) :: ngrid
55      character (len=*),intent(in) :: nom,titre,unite
56      integer,intent(in) :: dim
57      real,intent(in) :: px(ngrid,nbp_lev)
58
59! Local variables:
60
61      real*4 dx3(nbp_lon+1,nbp_lat,nbp_lev) ! to store a 3D data set
62      real*4 dx2(nbp_lon+1,nbp_lat)     ! to store a 2D (surface) data set
63      real*4 dx1(nbp_lev)           ! to store a 1D (column) data set
64      real*4 dx0
65      real*4 dx3_1d(1,nbp_lev) ! to store a profile with 1D model
66      real*4 dx2_1d ! to store a surface value with 1D model
67
68      real*4,save :: date
69!$OMP THREADPRIVATE(date)
70
71      REAL phis((nbp_lon+1),nbp_lat)
72      REAL area((nbp_lon+1),nbp_lat)
73
74      integer irythme
75      integer ierr,ierr2
76      integer i,j,l, 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,nbp_lev)
106      real dx3_glo(nbp_lon,nbp_lat,nbp_lev) ! to store a global 3D data set
107      real dx2_glop(klon_glo)
108      real dx2_glo(nbp_lon,nbp_lat)     ! to store a global 2D (surface) data set
109      real px2(ngrid)
110!      real dx1_glo(nbp_lev)          ! to store a 1D (column) data set
111!      real dx0_glo
112      real phisfi_glo(klon_glo) ! surface geopotential on global physics grid
113      real areafi_glo(klon_glo) ! mesh area on global physics grid
114#else
115      real phisfi_glo(ngrid) ! surface geopotential on global physics grid
116      real areafi_glo(ngrid) ! mesh area on global physics grid
117#endif
118
119!***************************************************************
120!Sortie des variables au rythme voulu
121
122      irythme = int(ecritphy) ! output rate set by ecritphy
123
124!***************************************************************
125
126! At very first call, check if there is a "diagfi.def" to use and read it
127! -----------------------------------------------------------------------
128      IF (firstcall) THEN
129         firstcall=.false.
130
131!$OMP MASTER
132  !      Open diagfi.def definition file if there is one:
133         open(99,file="diagfi.def",status='old',form='formatted',
134     s   iostat=ierr2)
135
136         if(ierr2.eq.0) then
137            diagfi_def=.true.
138            write(*,*) "******************"
139            write(*,*) "Reading diagfi.def"
140            write(*,*) "******************"
141            do n=1,n_nom_def_max
142              read(99,fmt='(a)',end=88) nom_def(n)
143              write(*,*) 'Output in diagfi: ', trim(nom_def(n))
144            end do
145 88         continue
146            if (n.ge.n_nom_def_max) then
147               write(*,*)"n_nom_def_max too small in writediagfi.F:",n
148               stop
149            end if
150            n_nom_def=n-1
151            close(99)
152         else
153            diagfi_def=.false.
154         endif
155!$OMP END MASTER
156!$OMP BARRIER
157      END IF ! of IF (firstcall)
158
159! Get out of write_diagfi if there is diagfi.def AND variable not listed
160!  ---------------------------------------------------------------------
161      if (diagfi_def) then
162          getout=.true.
163          do n=1,n_nom_def
164             if(trim(nom_def(n)).eq.nom) getout=.false.
165          end do
166          if (getout) return
167      end if
168
169! Initialisation of 'firstnom' and create/open the "diagfi.nc" NetCDF file
170! ------------------------------------------------------------------------
171! (at very first call to the subroutine, in accordance with diagfi.def)
172
173      if (firstnom.eq.'1234567890') then ! .true. for the very first valid
174      !   call to this subroutine; now set 'firstnom'
175         firstnom = nom
176         ! just to be sure, check that firstnom is large enough to hold nom
177         if (len_trim(firstnom).lt.len_trim(nom)) then
178           write(*,*) "writediagfi: Error !!!"
179           write(*,*) "   firstnom string not long enough!!"
180           write(*,*) "   increase its size to at least ",len_trim(nom)
181           stop
182         endif
183         
184         zitau = -1 ! initialize zitau
185
186#ifdef CPP_PARA
187          ! Gather phisfi() geopotential on physics grid
188          call Gather(phisfi,phisfi_glo)
189          ! Gather cell_area() mesh area on physics grid
190          call Gather(cell_area,areafi_glo)
191#else
192         phisfi_glo(:)=phisfi(:)
193         areafi_glo(:)=cell_area(:)
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         ! Build phis() and area()
221         IF (klon_glo>1) THEN
222          do i=1,nbp_lon+1 ! poles
223           phis(i,1)=phisfi_glo(1)
224           phis(i,nbp_lat)=phisfi_glo(klon_glo)
225           ! for area, divide at the poles by nbp_lon
226           area(i,1)=areafi_glo(1)/nbp_lon
227           area(i,nbp_lat)=areafi_glo(klon_glo)/nbp_lon
228          enddo
229          do j=2,nbp_lat-1
230           ig0= 1+(j-2)*nbp_lon
231           do i=1,nbp_lon
232              phis(i,j)=phisfi_glo(ig0+i)
233              area(i,j)=areafi_glo(ig0+i)
234           enddo
235           ! handle redundant point in longitude
236           phis(nbp_lon+1,j)=phis(1,j)
237           area(nbp_lon+1,j)=area(1,j)
238          enddo
239         ENDIF
240         
241         ! write "header" of file (longitudes, latitudes, geopotential, ...)
242         IF (klon_glo>1) THEN ! general 3D case
243           call iniwrite(nid,day_ini,phis,area,nbp_lon+1,nbp_lat)
244         ELSE ! 1D model
245           call iniwrite(nid,day_ini,phisfi_glo(1),areafi_glo(1),1,1)
246         ENDIF
247
248         endif ! of if (is_master)
249
250      else
251
252         if (is_master) then
253           ! only the master is required to do this
254
255           ! Open the NetCDF file
256           ierr = NF_OPEN(fichnom,NF_WRITE,nid)
257         endif ! of if (is_master)
258
259      endif ! if (firstnom.eq.'1234567890')
260
261! Increment time index 'zitau' if it is the "fist call" (at given time level)
262! to writediagfi
263!------------------------------------------------------------------------
264      if (nom.eq.firstnom) then
265          zitau = zitau + iphysiq
266      end if
267
268!--------------------------------------------------------
269! Write the variables to output file if it's time to do so
270!--------------------------------------------------------
271
272      if ( MOD(zitau+1,irythme) .eq.0.) then
273
274! Compute/write/extend 'Time' coordinate (date given in days)
275! (done every "first call" (at given time level) to writediagfi)
276! Note: date is incremented as 1 step ahead of physics time
277!--------------------------------------------------------
278
279        if (is_master) then
280           ! only the master is required to do this
281        if (nom.eq.firstnom) then
282        ! We have identified a "first call" (at given date)
283           ntime=ntime+1 ! increment # of stored time steps
284           ! compute corresponding date (in days and fractions thereof)
285           date= float (zitau +1)/float (day_step)
286           ! Get NetCDF ID of 'Time' variable
287           ierr= NF_INQ_VARID(nid,"Time",varid)
288           ! Write (append) the new date to the 'Time' array
289!#ifdef NC_DOUBLE
290!           ierr= NF_PUT_VARA_DOUBLE(nid,varid,ntime,1,date)
291!#else
292           ierr= NF_PUT_VARA_REAL(nid,varid,ntime,1,date)
293!#endif
294           if (ierr.ne.NF_NOERR) then
295              write(*,*) "***** PUT_VAR matter in writediagfi_nc"
296              write(*,*) "***** with time"
297              write(*,*) 'ierr=', ierr   
298c             call abort
299           endif
300
301           write(6,*)'WRITEDIAGFI: date= ', date
302        end if ! of if (nom.eq.firstnom)
303
304        endif ! of if (is_master)
305
306!Case of a 3D variable
307!---------------------
308        if (dim.eq.3) then
309
310#ifdef CPP_PARA
311          ! Gather field on a "global" (without redundant longitude) array
312          call Gather(px,dx3_glop)
313!$OMP MASTER
314          if (is_mpi_root) then
315            call Grid1Dto2D_glo(dx3_glop,dx3_glo)
316            ! copy dx3_glo() to dx3(:) and add redundant longitude
317            dx3(1:nbp_lon,:,:)=dx3_glo(1:nbp_lon,:,:)
318            dx3(nbp_lon+1,:,:)=dx3(1,:,:)
319          endif
320!$OMP END MASTER
321!$OMP BARRIER
322#else
323!         Passage variable physique -->  variable dynamique
324!         recast (copy) variable from physics grid to dynamics grid
325          IF (klon_glo>1) THEN ! General case
326           DO l=1,nbp_lev
327             DO i=1,nbp_lon+1
328                dx3(i,1,l)=px(1,l)
329                dx3(i,nbp_lat,l)=px(ngrid,l)
330             ENDDO
331             DO j=2,nbp_lat-1
332                ig0= 1+(j-2)*nbp_lon
333                DO i=1,nbp_lon
334                   dx3(i,j,l)=px(ig0+i,l)
335                ENDDO
336                dx3(nbp_lon+1,j,l)=dx3(1,j,l)
337             ENDDO
338           ENDDO
339          ELSE ! 1D model case
340           dx3_1d(1,1:nbp_lev)=px(1,1:nbp_lev)
341          ENDIF
342#endif
343!         Ecriture du champs
344
345          if (is_master) then
346           ! only the master writes to output
347! name of the variable
348           ierr= NF_INQ_VARID(nid,nom,varid)
349           if (ierr /= NF_NOERR) then
350! corresponding dimensions
351              ierr= NF_INQ_DIMID(nid,"longitude",id(1))
352              ierr= NF_INQ_DIMID(nid,"latitude",id(2))
353              ierr= NF_INQ_DIMID(nid,"altitude",id(3))
354              ierr= NF_INQ_DIMID(nid,"Time",id(4))
355
356! Create the variable if it doesn't exist yet
357
358              write (*,*) "=========================="
359              write (*,*) "DIAGFI: creating variable ",nom
360              call def_var(nid,nom,titre,unite,4,id,varid,ierr)
361
362           endif
363
364           corner(1)=1
365           corner(2)=1
366           corner(3)=1
367           corner(4)=ntime
368
369           IF (klon_glo==1) THEN
370             edges(1)=1
371           ELSE
372             edges(1)=nbp_lon+1
373           ENDIF
374           edges(2)=nbp_lat
375           edges(3)=nbp_lev
376           edges(4)=1
377!#ifdef NC_DOUBLE
378!           ierr= NF_PUT_VARA_DOUBLE(nid,varid,corner,edges,dx3)
379!#else
380!           write(*,*)"test:  nid=",nid," varid=",varid
381!           write(*,*)"       corner()=",corner
382!           write(*,*)"       edges()=",edges
383!           write(*,*)"       dx3()=",dx3
384           IF (klon_glo>1) THEN ! General case
385             ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx3)
386           ELSE
387             ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx3_1d)
388           ENDIF
389!#endif
390
391           if (ierr.ne.NF_NOERR) then
392              write(*,*) "***** PUT_VAR problem in writediagfi"
393              write(*,*) "***** with dx3: ",nom
394              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
395              stop
396           endif
397
398          endif !of if (is_master)
399
400!Case of a 2D variable
401!---------------------
402
403        else if (dim.eq.2) then
404
405#ifdef CPP_PARA
406          ! Gather field on a "global" (without redundant longitude) array
407          px2(:)=px(:,1)
408          call Gather(px2,dx2_glop)
409!$OMP MASTER
410          if (is_mpi_root) then
411            call Grid1Dto2D_glo(dx2_glop,dx2_glo)
412            ! copy dx2_glo() to dx2(:) and add redundant longitude
413            dx2(1:nbp_lon,:)=dx2_glo(1:nbp_lon,:)
414            dx2(nbp_lon+1,:)=dx2(1,:)
415          endif
416!$OMP END MASTER
417!$OMP BARRIER
418#else
419
420!         Passage variable physique -->  physique dynamique
421!         recast (copy) variable from physics grid to dynamics grid
422          IF (klon_glo>1) THEN ! General case
423             DO i=1,nbp_lon+1
424                dx2(i,1)=px(1,1)
425                dx2(i,nbp_lat)=px(ngrid,1)
426             ENDDO
427             DO j=2,nbp_lat-1
428                ig0= 1+(j-2)*nbp_lon
429                DO i=1,nbp_lon
430                   dx2(i,j)=px(ig0+i,1)
431                ENDDO
432                dx2(nbp_lon+1,j)=dx2(1,j)
433             ENDDO
434          ELSE ! 1D model case
435            dx2_1d=px(1,1)
436          ENDIF
437#endif
438
439          if (is_master) then
440           ! only the master writes to output
441!         write (*,*) 'In  writediagfi, on sauve:  ' , nom
442!         write (*,*) 'In  writediagfi. Estimated date = ' ,date
443           ierr= NF_INQ_VARID(nid,nom,varid)
444           if (ierr /= NF_NOERR) then
445! corresponding dimensions
446              ierr= NF_INQ_DIMID(nid,"longitude",id(1))
447              ierr= NF_INQ_DIMID(nid,"latitude",id(2))
448              ierr= NF_INQ_DIMID(nid,"Time",id(3))
449
450! Create the variable if it doesn't exist yet
451
452              write (*,*) "=========================="
453              write (*,*) "DIAGFI: creating variable ",nom
454
455              call def_var(nid,nom,titre,unite,3,id,varid,ierr)
456
457           endif
458
459           corner(1)=1
460           corner(2)=1
461           corner(3)=ntime
462           IF (klon_glo==1) THEN
463             edges(1)=1
464           ELSE
465             edges(1)=nbp_lon+1
466           ENDIF
467           edges(2)=nbp_lat
468           edges(3)=1
469
470
471!#ifdef NC_DOUBLE
472!           ierr = NF_PUT_VARA_DOUBLE (nid,varid,corner,edges,dx2)
473!#else         
474           IF (klon_glo>1) THEN ! General case
475             ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx2)
476           ELSE
477             ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx2_1d)
478           ENDIF
479!#endif     
480
481           if (ierr.ne.NF_NOERR) then
482              write(*,*) "***** PUT_VAR matter in writediagfi"
483              write(*,*) "***** with dx2: ",nom
484              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
485              stop
486           endif
487
488          endif !of if (is_master)
489
490!Case of a 1D variable (ie: a column)
491!---------------------------------------------------
492
493       else if (dim.eq.1) then
494         if (is_parallel) then
495           write(*,*) "writediagfi error: dim=1 not implemented ",
496     &                 "in parallel mode"
497           stop
498         endif
499!         Passage variable physique -->  physique dynamique
500!         recast (copy) variable from physics grid to dynamics grid
501          do l=1,nbp_lev
502            dx1(l)=px(1,l)
503          enddo
504         
505          ierr= NF_INQ_VARID(nid,nom,varid)
506           if (ierr /= NF_NOERR) then
507! corresponding dimensions
508              ierr= NF_INQ_DIMID(nid,"altitude",id(1))
509              ierr= NF_INQ_DIMID(nid,"Time",id(2))
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,2,id,varid,ierr)
517             
518           endif
519           
520           corner(1)=1
521           corner(2)=ntime
522           
523           edges(1)=nbp_lev
524           edges(2)=1
525!#ifdef NC_DOUBLE
526!           ierr= NF_PUT_VARA_DOUBLE(nid,varid,corner,edges,dx1)
527!#else
528           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx1)
529!#endif
530
531           if (ierr.ne.NF_NOERR) then
532              write(*,*) "***** PUT_VAR problem in writediagfi"
533              write(*,*) "***** with dx1: ",nom
534              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
535              stop
536           endif
537
538!Case of a 0D variable (ie: a time-dependent scalar)
539!---------------------------------------------------
540
541        else if (dim.eq.0) then
542
543           dx0 = px (1,1)
544
545          if (is_master) then
546           ! only the master writes to output
547           ierr= NF_INQ_VARID(nid,nom,varid)
548           if (ierr /= NF_NOERR) then
549! corresponding dimensions
550              ierr= NF_INQ_DIMID(nid,"Time",id(1))
551
552! Create the variable if it doesn't exist yet
553
554              write (*,*) "=========================="
555              write (*,*) "DIAGFI: creating variable ",nom
556
557              call def_var(nid,nom,titre,unite,1,id,varid,ierr)
558
559           endif
560
561           corner(1)=ntime
562           edges(1)=1
563
564!#ifdef NC_DOUBLE
565!           ierr = NF_PUT_VARA_DOUBLE (nid,varid,corner,edges,dx0) 
566!#else
567           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx0)
568!#endif
569           if (ierr.ne.NF_NOERR) then
570              write(*,*) "***** PUT_VAR matter in writediagfi"
571              write(*,*) "***** with dx0: ",nom
572              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
573              stop
574           endif
575
576          endif !of if (is_master)
577
578        endif ! of if (dim.eq.3) elseif(dim.eq.2)...
579
580      endif ! of if ( MOD(zitau+1,irythme) .eq.0.)
581
582      if (is_master) then
583        ierr= NF_CLOSE(nid)
584      endif
585
586#endif
587! of #ifndef MESOSCALE
588      end
Note: See TracBrowser for help on using the repository browser.