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

Last change on this file since 2443 was 2311, checked in by emillour, 5 years ago

Mars GCM:
Code tidying: use getin_p() instead of getin() and use "call abort_physic"
instead of "stop" or "call abort".
EM

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