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

Last change on this file since 1150 was 993, checked in by emillour, 11 years ago

Generic GCM:

  • Some more cleanup in dynamics:
    • Moved "start2archive" (and auxilliary routines) to phystd
    • removed unused (obsolete) testharm.F , para_netcdf.h , readhead_NC.F , angtot.h from dyn3d
    • removed obsolete addit.F (and change corresponding lines in gcm)
    • remove unused "description.h" (and many places where it was "included")

EM

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