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

Last change on this file since 1525 was 1525, checked in by emillour, 9 years ago

All GCMs:
More on enforcing dynamics/physics separation: get rid of references to "control_mod" from physics packages.
EM

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