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

Last change on this file since 1524 was 1524, checked in by emillour, 9 years ago

All GCMS:
More updates to enforce dynamics/physics separation:

get rid of references to "temps_mod" from physics packages;
make a "time_phylmdz_mod.F90" module to store that
information and fill it via "iniphysiq".

EM

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