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

Last change on this file since 2156 was 2141, checked in by emillour, 5 years ago

Mars GCM:

  • minor updates: enable output of Ls in diagfi files and make an error message more verbose.

EM

File size: 20.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#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=(zitau +1.)/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,": ",NF_STRERROR(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 ",trim(nom)
358              call def_var(nid,nom,titre,unite,4,id,varid,ierr)
359
360           else
361             if (ntime==0) then
362              write(*,*) "DIAGFI Error: failed creating variable ",
363     &                   trim(nom)
364              write(*,*) "it seems it already exists!"
365              stop
366             endif
367           endif
368
369           corner(1)=1
370           corner(2)=1
371           corner(3)=1
372           corner(4)=ntime
373
374           IF (klon_glo==1) THEN
375             edges(1)=1
376           ELSE
377             edges(1)=nbp_lon+1
378           ENDIF
379           edges(2)=nbp_lat
380           edges(3)=nbp_lev
381           edges(4)=1
382!#ifdef NC_DOUBLE
383!           ierr= NF_PUT_VARA_DOUBLE(nid,varid,corner,edges,dx3)
384!#else
385!           write(*,*)"test:  nid=",nid," varid=",varid
386!           write(*,*)"       corner()=",corner
387!           write(*,*)"       edges()=",edges
388!           write(*,*)"       dx3()=",dx3
389           IF (klon_glo>1) THEN ! General case
390             ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx3)
391           ELSE
392             ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx3_1d)
393           ENDIF
394!#endif
395
396           if (ierr.ne.NF_NOERR) then
397              write(*,*) "***** PUT_VAR problem in writediagfi"
398              write(*,*) "***** with dx3: ",trim(nom)
399              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
400              stop
401           endif
402
403          endif !of if (is_master)
404
405!Case of a 2D variable
406!---------------------
407
408        else if (dim.eq.2) then
409
410#ifdef CPP_PARA
411          ! Gather field on a "global" (without redundant longitude) array
412          px2(:)=px(:,1)
413          call Gather(px2,dx2_glop)
414!$OMP MASTER
415          if (is_mpi_root) then
416            call Grid1Dto2D_glo(dx2_glop,dx2_glo)
417            ! copy dx2_glo() to dx2(:) and add redundant longitude
418            dx2(1:nbp_lon,:)=dx2_glo(1:nbp_lon,:)
419            dx2(nbp_lon+1,:)=dx2(1,:)
420          endif
421!$OMP END MASTER
422!$OMP BARRIER
423#else
424
425!         Passage variable physique -->  physique dynamique
426!         recast (copy) variable from physics grid to dynamics grid
427          IF (klon_glo>1) THEN ! General case
428             DO i=1,nbp_lon+1
429                dx2(i,1)=px(1,1)
430                dx2(i,nbp_lat)=px(ngrid,1)
431             ENDDO
432             DO j=2,nbp_lat-1
433                ig0= 1+(j-2)*nbp_lon
434                DO i=1,nbp_lon
435                   dx2(i,j)=px(ig0+i,1)
436                ENDDO
437                dx2(nbp_lon+1,j)=dx2(1,j)
438             ENDDO
439          ELSE ! 1D model case
440            dx2_1d=px(1,1)
441          ENDIF
442#endif
443
444          if (is_master) then
445           ! only the master writes to output
446!         write (*,*) 'In  writediagfi, on sauve:  ' , nom
447!         write (*,*) 'In  writediagfi. Estimated date = ' ,date
448           ierr= NF_INQ_VARID(nid,nom,varid)
449           if (ierr /= NF_NOERR) then
450! corresponding dimensions
451              ierr= NF_INQ_DIMID(nid,"longitude",id(1))
452              ierr= NF_INQ_DIMID(nid,"latitude",id(2))
453              ierr= NF_INQ_DIMID(nid,"Time",id(3))
454
455! Create the variable if it doesn't exist yet
456
457              write (*,*) "=========================="
458              write (*,*) "DIAGFI: creating variable ",trim(nom)
459
460              call def_var(nid,nom,titre,unite,3,id,varid,ierr)
461
462           else
463             if (ntime==0) then
464              write(*,*) "DIAGFI Error: failed creating variable ",
465     &                   trim(nom)
466              write(*,*) "it seems it already exists!"
467              stop
468             endif
469           endif
470
471           corner(1)=1
472           corner(2)=1
473           corner(3)=ntime
474           IF (klon_glo==1) THEN
475             edges(1)=1
476           ELSE
477             edges(1)=nbp_lon+1
478           ENDIF
479           edges(2)=nbp_lat
480           edges(3)=1
481
482
483!#ifdef NC_DOUBLE
484!           ierr = NF_PUT_VARA_DOUBLE (nid,varid,corner,edges,dx2)
485!#else         
486           IF (klon_glo>1) THEN ! General case
487             ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx2)
488           ELSE
489             ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx2_1d)
490           ENDIF
491!#endif     
492
493           if (ierr.ne.NF_NOERR) then
494              write(*,*) "***** PUT_VAR matter in writediagfi"
495              write(*,*) "***** with dx2: ",trim(nom)
496              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
497              stop
498           endif
499
500          endif !of if (is_master)
501
502!Case of a 1D variable (ie: a column)
503!---------------------------------------------------
504
505       else if (dim.eq.1) then
506         if (is_parallel) then
507           write(*,*) "writediagfi error: dim=1 not implemented ",
508     &                 "in parallel mode. Problem for ",trim(nom)
509           stop
510         endif
511!         Passage variable physique -->  physique dynamique
512!         recast (copy) variable from physics grid to dynamics grid
513          do l=1,nbp_lev
514            dx1(l)=px(1,l)
515          enddo
516         
517          ierr= NF_INQ_VARID(nid,nom,varid)
518           if (ierr /= NF_NOERR) then
519! corresponding dimensions
520              ierr= NF_INQ_DIMID(nid,"altitude",id(1))
521              ierr= NF_INQ_DIMID(nid,"Time",id(2))
522
523! Create the variable if it doesn't exist yet
524
525              write (*,*) "=========================="
526              write (*,*) "DIAGFI: creating variable ",trim(nom)
527
528              call def_var(nid,nom,titre,unite,2,id,varid,ierr)
529             
530           else
531             if (ntime==0) then
532              write(*,*) "DIAGFI Error: failed creating variable ",
533     &                   trim(nom)
534              write(*,*) "it seems it already exists!"
535              stop
536             endif
537           endif
538           
539           corner(1)=1
540           corner(2)=ntime
541           
542           edges(1)=nbp_lev
543           edges(2)=1
544!#ifdef NC_DOUBLE
545!           ierr= NF_PUT_VARA_DOUBLE(nid,varid,corner,edges,dx1)
546!#else
547           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx1)
548!#endif
549
550           if (ierr.ne.NF_NOERR) then
551              write(*,*) "***** PUT_VAR problem in writediagfi"
552              write(*,*) "***** with dx1: ",trim(nom)
553              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
554              stop
555           endif
556
557!Case of a 0D variable (ie: a time-dependent scalar)
558!---------------------------------------------------
559
560        else if (dim.eq.0) then
561
562           dx0 = px (1,1)
563
564          if (is_master) then
565           ! only the master writes to output
566           ierr= NF_INQ_VARID(nid,nom,varid)
567           if (ierr /= NF_NOERR) then
568! corresponding dimensions
569              ierr= NF_INQ_DIMID(nid,"Time",id(1))
570
571! Create the variable if it doesn't exist yet
572
573              write (*,*) "=========================="
574              write (*,*) "DIAGFI: creating variable ",trim(nom)
575
576              call def_var(nid,nom,titre,unite,1,id,varid,ierr)
577
578           else
579             if (ntime==0) then
580              write(*,*) "DIAGFI Error: failed creating variable ",
581     &                   trim(nom)
582              write(*,*) "it seems it already exists!"
583              stop
584             endif
585           endif
586
587           corner(1)=ntime
588           edges(1)=1
589
590!#ifdef NC_DOUBLE
591!           ierr = NF_PUT_VARA_DOUBLE (nid,varid,corner,edges,dx0) 
592!#else
593           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx0)
594!#endif
595           if (ierr.ne.NF_NOERR) then
596              write(*,*) "***** PUT_VAR matter in writediagfi"
597              write(*,*) "***** with dx0: ",trim(nom)
598              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
599              stop
600           endif
601
602          endif !of if (is_master)
603
604        endif ! of if (dim.eq.3) elseif(dim.eq.2)...
605
606      endif ! of if ( MOD(zitau+1,irythme) .eq.0.)
607
608      if (is_master) then
609        ierr= NF_CLOSE(nid)
610      endif
611
612      end
Note: See TracBrowser for help on using the repository browser.