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

Last change on this file since 1421 was 1318, checked in by aslmd, 11 years ago

LMDZ.GENERIC. Two fixes for 1D version. First a conflict between comvert.h and planete_mod.F for ap bp preff (this one should be given a closer look, this redundancy is weird). Secondly output frequency was hardcoded to 1 (instead of ecritphy) in r1216.

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