source: trunk/LMDZ.MARS/libf/phymars/writediagfi.F @ 1532

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

Mars GCM:

  • Some fixes for buggy outputs in 1D introduced by previous code modifications.
  • made statto.h a module.
  • ecritphy in dyn3d/control_mod.F90 should be an integer, not a real.

EM

File size: 19.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 comgeomphy, only: airephy
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#ifdef CPP_PARA
102! Added to work in parallel mode
103      real dx3_glop(klon_glo,nbp_lev)
104      real dx3_glo(nbp_lon,nbp_lat,nbp_lev) ! to store a global 3D data set
105      real dx2_glop(klon_glo)
106      real dx2_glo(nbp_lon,nbp_lat)     ! to store a global 2D (surface) data set
107      real px2(ngrid)
108!      real dx1_glo(nbp_lev)          ! to store a 1D (column) data set
109!      real dx0_glo
110      real phisfi_glo(klon_glo) ! surface geopotential on global physics grid
111      real areafi_glo(klon_glo) ! mesh area on global physics grid
112#else
113      real phisfi_glo(ngrid) ! surface geopotential on global physics grid
114      real areafi_glo(ngrid) ! mesh area on global physics grid
115#endif
116
117!***************************************************************
118!Sortie des variables au rythme voulu
119
120      irythme = int(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          ! Gather airephy() mesh area on physics grid
192          call Gather(airephy,areafi_glo)
193#else
194         phisfi_glo(:)=phisfi(:)
195         areafi_glo(:)=airephy(:)
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         ! Build phis() and area()
223         IF (klon_glo>1) THEN
224          do i=1,nbp_lon+1 ! poles
225           phis(i,1)=phisfi_glo(1)
226           phis(i,nbp_lat)=phisfi_glo(klon_glo)
227           ! for area, divide at the poles by nbp_lon
228           area(i,1)=areafi_glo(1)/nbp_lon
229           area(i,nbp_lat)=areafi_glo(klon_glo)/nbp_lon
230          enddo
231          do j=2,nbp_lat-1
232           ig0= 1+(j-2)*nbp_lon
233           do i=1,nbp_lon
234              phis(i,j)=phisfi_glo(ig0+i)
235              area(i,j)=areafi_glo(ig0+i)
236           enddo
237           ! handle redundant point in longitude
238           phis(nbp_lon+1,j)=phis(1,j)
239           area(nbp_lon+1,j)=area(1,j)
240          enddo
241         ENDIF
242         
243         ! write "header" of file (longitudes, latitudes, geopotential, ...)
244         IF (klon_glo>1) THEN ! general 3D case
245           call iniwrite(nid,day_ini,phis,area,nbp_lon+1,nbp_lat)
246         ELSE ! 1D model
247           call iniwrite(nid,day_ini,phisfi_glo(1),areafi_glo(1),1,1)
248         ENDIF
249
250         endif ! of if (is_master)
251
252      else
253
254         if (is_master) then
255           ! only the master is required to do this
256
257           ! Open the NetCDF file
258           ierr = NF_OPEN(fichnom,NF_WRITE,nid)
259         endif ! of if (is_master)
260
261      endif ! if (firstnom.eq.'1234567890')
262
263      if (klon_glo.eq.1) then
264        ! in testphys1d, for the 1d version of the GCM, iphysiq and irythme
265        ! are undefined; so set them to 1
266        iphysiq=1
267        irythme=1
268      endif
269
270! Increment time index 'zitau' if it is the "fist call" (at given time level)
271! to writediagfi
272!------------------------------------------------------------------------
273      if (nom.eq.firstnom) then
274          zitau = zitau + iphysiq
275      end if
276
277!--------------------------------------------------------
278! Write the variables to output file if it's time to do so
279!--------------------------------------------------------
280
281      if ( MOD(zitau+1,irythme) .eq.0.) then
282
283! Compute/write/extend 'Time' coordinate (date given in days)
284! (done every "first call" (at given time level) to writediagfi)
285! Note: date is incremented as 1 step ahead of physics time
286!--------------------------------------------------------
287
288        if (is_master) then
289           ! only the master is required to do this
290        if (nom.eq.firstnom) then
291        ! We have identified a "first call" (at given date)
292           ntime=ntime+1 ! increment # of stored time steps
293           ! compute corresponding date (in days and fractions thereof)
294           date= float (zitau +1)/float (day_step)
295           ! Get NetCDF ID of 'Time' variable
296           ierr= NF_INQ_VARID(nid,"Time",varid)
297           ! Write (append) the new date to the 'Time' array
298!#ifdef NC_DOUBLE
299!           ierr= NF_PUT_VARA_DOUBLE(nid,varid,ntime,1,date)
300!#else
301           ierr= NF_PUT_VARA_REAL(nid,varid,ntime,1,date)
302!#endif
303           if (ierr.ne.NF_NOERR) then
304              write(*,*) "***** PUT_VAR matter in writediagfi_nc"
305              write(*,*) "***** with time"
306              write(*,*) 'ierr=', ierr   
307c             call abort
308           endif
309
310           write(6,*)'WRITEDIAGFI: date= ', date
311        end if ! of if (nom.eq.firstnom)
312
313        endif ! of if (is_master)
314
315!Case of a 3D variable
316!---------------------
317        if (dim.eq.3) then
318
319#ifdef CPP_PARA
320          ! Gather field on a "global" (without redundant longitude) array
321          call Gather(px,dx3_glop)
322!$OMP MASTER
323          if (is_mpi_root) then
324            call Grid1Dto2D_glo(dx3_glop,dx3_glo)
325            ! copy dx3_glo() to dx3(:) and add redundant longitude
326            dx3(1:nbp_lon,:,:)=dx3_glo(1:nbp_lon,:,:)
327            dx3(nbp_lon+1,:,:)=dx3(1,:,:)
328          endif
329!$OMP END MASTER
330!$OMP BARRIER
331#else
332!         Passage variable physique -->  variable dynamique
333!         recast (copy) variable from physics grid to dynamics grid
334          IF (klon_glo>1) THEN ! General case
335           DO l=1,nbp_lev
336             DO i=1,nbp_lon+1
337                dx3(i,1,l)=px(1,l)
338                dx3(i,nbp_lat,l)=px(ngrid,l)
339             ENDDO
340             DO j=2,nbp_lat-1
341                ig0= 1+(j-2)*nbp_lon
342                DO i=1,nbp_lon
343                   dx3(i,j,l)=px(ig0+i,l)
344                ENDDO
345                dx3(nbp_lon+1,j,l)=dx3(1,j,l)
346             ENDDO
347           ENDDO
348          ELSE ! 1D model case
349           dx3_1d(1,1:nbp_lev)=px(1,1:nbp_lev)
350          ENDIF
351#endif
352!         Ecriture du champs
353
354          if (is_master) then
355           ! only the master writes to output
356! name of the variable
357           ierr= NF_INQ_VARID(nid,nom,varid)
358           if (ierr /= NF_NOERR) then
359! corresponding dimensions
360              ierr= NF_INQ_DIMID(nid,"longitude",id(1))
361              ierr= NF_INQ_DIMID(nid,"latitude",id(2))
362              ierr= NF_INQ_DIMID(nid,"altitude",id(3))
363              ierr= NF_INQ_DIMID(nid,"Time",id(4))
364
365! Create the variable if it doesn't exist yet
366
367              write (*,*) "=========================="
368              write (*,*) "DIAGFI: creating variable ",nom
369              call def_var(nid,nom,titre,unite,4,id,varid,ierr)
370
371           endif
372
373           corner(1)=1
374           corner(2)=1
375           corner(3)=1
376           corner(4)=ntime
377
378           IF (klon_glo==1) THEN
379             edges(1)=1
380           ELSE
381             edges(1)=nbp_lon+1
382           ENDIF
383           edges(2)=nbp_lat
384           edges(3)=nbp_lev
385           edges(4)=1
386!#ifdef NC_DOUBLE
387!           ierr= NF_PUT_VARA_DOUBLE(nid,varid,corner,edges,dx3)
388!#else
389!           write(*,*)"test:  nid=",nid," varid=",varid
390!           write(*,*)"       corner()=",corner
391!           write(*,*)"       edges()=",edges
392!           write(*,*)"       dx3()=",dx3
393           IF (klon_glo>1) THEN ! General case
394             ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx3)
395           ELSE
396             ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx3_1d)
397           ENDIF
398!#endif
399
400           if (ierr.ne.NF_NOERR) then
401              write(*,*) "***** PUT_VAR problem in writediagfi"
402              write(*,*) "***** with dx3: ",nom
403              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
404              stop
405           endif
406
407          endif !of if (is_master)
408
409!Case of a 2D variable
410!---------------------
411
412        else if (dim.eq.2) then
413
414#ifdef CPP_PARA
415          ! Gather field on a "global" (without redundant longitude) array
416          px2(:)=px(:,1)
417          call Gather(px2,dx2_glop)
418!$OMP MASTER
419          if (is_mpi_root) then
420            call Grid1Dto2D_glo(dx2_glop,dx2_glo)
421            ! copy dx2_glo() to dx2(:) and add redundant longitude
422            dx2(1:nbp_lon,:)=dx2_glo(1:nbp_lon,:)
423            dx2(nbp_lon+1,:)=dx2(1,:)
424          endif
425!$OMP END MASTER
426!$OMP BARRIER
427#else
428
429!         Passage variable physique -->  physique dynamique
430!         recast (copy) variable from physics grid to dynamics grid
431          IF (klon_glo>1) THEN ! General case
432             DO i=1,nbp_lon+1
433                dx2(i,1)=px(1,1)
434                dx2(i,nbp_lat)=px(ngrid,1)
435             ENDDO
436             DO j=2,nbp_lat-1
437                ig0= 1+(j-2)*nbp_lon
438                DO i=1,nbp_lon
439                   dx2(i,j)=px(ig0+i,1)
440                ENDDO
441                dx2(nbp_lon+1,j)=dx2(1,j)
442             ENDDO
443          ELSE ! 1D model case
444            dx2_1d=px(1,1)
445          ENDIF
446#endif
447
448          if (is_master) then
449           ! only the master writes to output
450!         write (*,*) 'In  writediagfi, on sauve:  ' , nom
451!         write (*,*) 'In  writediagfi. Estimated date = ' ,date
452           ierr= NF_INQ_VARID(nid,nom,varid)
453           if (ierr /= NF_NOERR) then
454! corresponding dimensions
455              ierr= NF_INQ_DIMID(nid,"longitude",id(1))
456              ierr= NF_INQ_DIMID(nid,"latitude",id(2))
457              ierr= NF_INQ_DIMID(nid,"Time",id(3))
458
459! Create the variable if it doesn't exist yet
460
461              write (*,*) "=========================="
462              write (*,*) "DIAGFI: creating variable ",nom
463
464              call def_var(nid,nom,titre,unite,3,id,varid,ierr)
465
466           endif
467
468           corner(1)=1
469           corner(2)=1
470           corner(3)=ntime
471           IF (klon_glo==1) THEN
472             edges(1)=1
473           ELSE
474             edges(1)=nbp_lon+1
475           ENDIF
476           edges(2)=nbp_lat
477           edges(3)=1
478
479
480!#ifdef NC_DOUBLE
481!           ierr = NF_PUT_VARA_DOUBLE (nid,varid,corner,edges,dx2)
482!#else         
483           IF (klon_glo>1) THEN ! General case
484             ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx2)
485           ELSE
486             ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx2_1d)
487           ENDIF
488!#endif     
489
490           if (ierr.ne.NF_NOERR) then
491              write(*,*) "***** PUT_VAR matter in writediagfi"
492              write(*,*) "***** with dx2: ",nom
493              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
494              stop
495           endif
496
497          endif !of if (is_master)
498
499!Case of a 1D variable (ie: a column)
500!---------------------------------------------------
501
502       else if (dim.eq.1) then
503         if (is_parallel) then
504           write(*,*) "writediagfi error: dim=1 not implemented ",
505     &                 "in parallel mode"
506           stop
507         endif
508!         Passage variable physique -->  physique dynamique
509!         recast (copy) variable from physics grid to dynamics grid
510          do l=1,nbp_lev
511            dx1(l)=px(1,l)
512          enddo
513         
514          ierr= NF_INQ_VARID(nid,nom,varid)
515           if (ierr /= NF_NOERR) then
516! corresponding dimensions
517              ierr= NF_INQ_DIMID(nid,"altitude",id(1))
518              ierr= NF_INQ_DIMID(nid,"Time",id(2))
519
520! Create the variable if it doesn't exist yet
521
522              write (*,*) "=========================="
523              write (*,*) "DIAGFI: creating variable ",nom
524
525              call def_var(nid,nom,titre,unite,2,id,varid,ierr)
526             
527           endif
528           
529           corner(1)=1
530           corner(2)=ntime
531           
532           edges(1)=nbp_lev
533           edges(2)=1
534!#ifdef NC_DOUBLE
535!           ierr= NF_PUT_VARA_DOUBLE(nid,varid,corner,edges,dx1)
536!#else
537           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx1)
538!#endif
539
540           if (ierr.ne.NF_NOERR) then
541              write(*,*) "***** PUT_VAR problem in writediagfi"
542              write(*,*) "***** with dx1: ",nom
543              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
544              stop
545           endif
546
547!Case of a 0D variable (ie: a time-dependent scalar)
548!---------------------------------------------------
549
550        else if (dim.eq.0) then
551
552           dx0 = px (1,1)
553
554          if (is_master) then
555           ! only the master writes to output
556           ierr= NF_INQ_VARID(nid,nom,varid)
557           if (ierr /= NF_NOERR) then
558! corresponding dimensions
559              ierr= NF_INQ_DIMID(nid,"Time",id(1))
560
561! Create the variable if it doesn't exist yet
562
563              write (*,*) "=========================="
564              write (*,*) "DIAGFI: creating variable ",nom
565
566              call def_var(nid,nom,titre,unite,1,id,varid,ierr)
567
568           endif
569
570           corner(1)=ntime
571           edges(1)=1
572
573!#ifdef NC_DOUBLE
574!           ierr = NF_PUT_VARA_DOUBLE (nid,varid,corner,edges,dx0) 
575!#else
576           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx0)
577!#endif
578           if (ierr.ne.NF_NOERR) then
579              write(*,*) "***** PUT_VAR matter in writediagfi"
580              write(*,*) "***** with dx0: ",nom
581              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
582              stop
583           endif
584
585          endif !of if (is_master)
586
587        endif ! of if (dim.eq.3) elseif(dim.eq.2)...
588
589      endif ! of if ( MOD(zitau+1,irythme) .eq.0.)
590
591      if (is_master) then
592        ierr= NF_CLOSE(nid)
593      endif
594
595      end
Note: See TracBrowser for help on using the repository browser.