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

Last change on this file since 1151 was 993, checked in by emillour, 12 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
RevLine 
[135]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)
[323]8!  (ou encore 1d, dans le cas de testphys1d, pour sortir une colonne)
[135]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
[323]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
[135]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!=================================================================
[787]41
42      USE surfdat_h
[965]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
[135]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:
[323]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)
[135]65
66! Local variables:
67
[323]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
[135]72
[323]73      real*4,save :: date
[135]74
[965]75      REAL phis(ip1jmp1) ! surface geopotential on extended lonxlat grid
[135]76
77      integer irythme
[323]78      integer ierr,ierr2
[135]79      integer iq
80      integer i,j,l,zmax , ig0
81
[323]82      integer,save :: zitau=0
83      character(len=20),save :: firstnom='1234567890'
[135]84
85! Ajouts
86      integer, save :: ntime=0
87      integer :: idim,varid
88      integer :: nid
[323]89      character(len=*),parameter :: fichnom="diagfi.nc"
[135]90      integer, dimension(4) :: id
91      integer, dimension(4) :: edges,corner
[323]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.
[787]101
[862]102      integer dimvert
103
[965]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
[135]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
[787]131
[323]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.
[135]136
[787]137!      Open diagfi.def definition file if there is one:
[323]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
[862]172
[135]173! Initialisation of 'firstnom' and create/open the "diagfi.nc" NetCDF file
174! ------------------------------------------------------------------------
[323]175! (at very first call to the subroutine, in accordance with diagfi.def)
[135]176
[323]177      if (firstnom.eq.'1234567890') then ! .true. for the very first valid
178      !   call to this subroutine; now set 'firstnom'
[135]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
[862]188         zitau = -1 ! initialize zitau
189
190#ifdef CPP_PARA
[965]191          ! Gather phisfi() geopotential on physics grid
192          call Gather(phisfi,phisfi_glo)
193#else
194         phisfi_glo(:)=phisfi(:)
195#endif
196
[862]197         !! parallel: we cannot use the usual writediagfi method
[965]198!!         call iophys_ini
199         if (is_master) then
200         ! only the master is required to do this
201
[135]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
[323]207!#ifdef NC_DOUBLE
208!         ierr = NF_DEF_VAR (nid, "Time", NF_DOUBLE, 1, idim,varid)
209!#else
[135]210         ierr = NF_DEF_VAR (nid, "Time", NF_FLOAT, 1, idim,varid)
[323]211!#endif
[135]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, ...)
[965]222         call gr_fi_dyn(1,size(phisfi_glo),iip1,jjp1,phisfi_glo,phis)
[135]223         call iniwrite(nid,day_ini,phis)
[862]224
[965]225         endif ! of if (is_master)
226
[135]227      else
[965]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
[135]236      endif ! if (firstnom.eq.'1234567890')
237
[787]238      if (ngrid.eq.1) then
[323]239        ! in testphys1d, for the 1d version of the GCM, iphysiq and irythme
[135]240        ! are undefined; so set them to 1
241        iphysiq=1
[526]242        !JL11 irythme=1
[135]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
[965]264        if (is_master) then
265           ! only the master is required to do this
[135]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
[323]274!#ifdef NC_DOUBLE
275!           ierr= NF_PUT_VARA_DOUBLE(nid,varid,ntime,1,date)
276!#else
[135]277           ierr= NF_PUT_VARA_REAL(nid,varid,ntime,1,date)
[323]278!#endif
[135]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
[965]289        endif ! of if (is_master)
290
[135]291!Case of a 3D variable
292!---------------------
293        if (dim.eq.3) then
294
[965]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
[135]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
[965]323#endif
[135]324!         Ecriture du champs
325
[965]326          if (is_master) then
327           ! only the master writes to output
[135]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
[323]354!#ifdef NC_DOUBLE
355!           ierr= NF_PUT_VARA_DOUBLE(nid,varid,corner,edges,dx3)
356!#else
[965]357!           write(*,*)"test:  nid=",nid," varid=",varid
358!           write(*,*)"       corner()=",corner
359!           write(*,*)"       edges()=",edges
360!           write(*,*)"       dx3()=",dx3
[135]361           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx3)
[323]362!#endif
[135]363
364           if (ierr.ne.NF_NOERR) then
365              write(*,*) "***** PUT_VAR problem in writediagfi"
366              write(*,*) "***** with ",nom
[965]367              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
[135]368c             call abort
369           endif
370
[965]371          endif !of if (is_master)
372
[135]373!Case of a 2D variable
374!---------------------
375
376        else if (dim.eq.2) then
377
[965]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
[135]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
[965]407#endif
[135]408
[965]409          if (is_master) then
410           ! only the master writes to output
[135]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
[323]437!#ifdef NC_DOUBLE
438!           ierr = NF_PUT_VARA_DOUBLE (nid,varid,corner,edges,dx2)
439!#else         
[135]440           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx2)
[323]441!#endif     
[135]442
443           if (ierr.ne.NF_NOERR) then
444              write(*,*) "***** PUT_VAR matter in writediagfi"
445              write(*,*) "***** with ",nom
[965]446              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
[135]447c             call abort
448           endif
449
[965]450          endif !of if (is_master)
451
[135]452!Case of a 1D variable (ie: a column)
453!---------------------------------------------------
454
455       else if (dim.eq.1) then
[965]456         if (is_parallel) then
457           write(*,*) "writediagfi error: dim=1 not implemented ",
458     &                 "in parallel mode"
459           stop
460         endif
[135]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
[323]487!#ifdef NC_DOUBLE
488!           ierr= NF_PUT_VARA_DOUBLE(nid,varid,corner,edges,dx1)
489!#else
[135]490           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx1)
[323]491!#endif
[135]492
493           if (ierr.ne.NF_NOERR) then
494              write(*,*) "***** PUT_VAR problem in writediagfi"
495              write(*,*) "***** with ",nom
[969]496              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
[135]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
[965]504
[135]505           dx0 = px (1,1)
506
[969]507          if (is_master) then
508           ! only the master writes to output
[135]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
[323]526!#ifdef NC_DOUBLE
527!           ierr = NF_PUT_VARA_DOUBLE (nid,varid,corner,edges,dx0) 
528!#else
[135]529           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx0)
[323]530!#endif
[135]531           if (ierr.ne.NF_NOERR) then
532              write(*,*) "***** PUT_VAR matter in writediagfi"
533              write(*,*) "***** with ",nom
[969]534              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
[135]535c             call abort
536           endif
537
[969]538          endif !of if (is_master)
539
[135]540        endif ! of if (dim.eq.3) elseif(dim.eq.2)...
541
542      endif ! of if ( MOD(zitau+1,irythme) .eq.0.)
543
[965]544      if (is_master) then
545        ierr= NF_CLOSE(nid)
546      endif
[787]547
[135]548      end
[965]549
Note: See TracBrowser for help on using the repository browser.