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

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

Mars GCM:

  • fix for 1D in writediagfi to enable writing at "ecritphy" rate.
  • removed iniprint.h from phymars/dyn1d since it is in "misc"
  • Some code cleanup in anticipation of future updates:
    • changed variable names in comgeomphy.F90: give them more explicit names: rlond => longitude , rlatd => latitude, airephy => cell_area, cuphy => dx , cvphy => dy
    • removed long(), lati() and area() from comgeomfi_h.F90, use longitude(), latitude() and cell_are() from comgeomphy.F90 instead

EM

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