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

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

Generic GCM: Towards a cleaner separation between physics and dynamics

  • Got rid of references to "dimensions.h" from physics packages: use nbp_lon (=iim), nbp_lat (=jjp1) and nbp_lev from module mod_grid_phy_lmdz (in phy_common) instead.
  • Removed module "comhdiff_mod.F90", as it is only used by module surf_heat_transp_mod.F90, moved module variables there.
  • Addedin "surf_heat_transp_mod" local versions of some arrays and routines (from dyn3d) required to compute gradient, divergence, etc. on the global dynamics grid. As before, the slab ocean only works in serial.

EM

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