source: trunk/LMDZ.GENERIC/libf/phystd/writediagspecVI.F @ 2733

Last change on this file since 2733 was 2732, checked in by aslmd, 2 years ago

took into account unstructured grid in writediagspec*.F

File size: 11.3 KB
Line 
1      subroutine writediagspecVI(ngrid,nom,titre,unite,dimpx,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!  Dans la version 2000, la periode d'ecriture est celle de
9!  "ecritphy " regle dans le fichier de controle de run :  run.def
10!
11!    writediagfi peut etre appele de n'importe quelle subroutine
12!    de la physique, plusieurs fois. L'initialisation et la creation du
13!    fichier se fait au tout premier appel.
14!
15! WARNING : les variables dynamique (u,v,t,q,ps)
16!  sauvees par writediagfi avec une
17! date donnee sont legerement differentes que dans le fichier histoire car
18! on ne leur a pas encore ajoute de la dissipation et de la physique !!!
19! IL est  RECOMMANDE d'ajouter les tendance physique a ces variables
20! avant l'ecriture dans diagfi (cf. physiq.F)
21
22!
23!  parametres (input) :
24!  ----------
25!      ngrid : nombres de point ou est calcule la physique
26!                (ngrid = 2+(jjm-1)*iim - 1/jjm)
27!                 (= nlon ou klon dans la physique terrestre)
28!     
29!      unit : unite logique du fichier de sortie (toujours la meme)
30!      nom  : nom de la variable a sortir (chaine de caracteres)
31!      titre: titre de la variable (chaine de caracteres)
32!      unite : unite de la variable (chaine de caracteres)
33!      px : variable a sortir (real 0, 2, ou 3d)
34!      dimpx : dimension de px : 0, 2, ou 3 dimensions
35!
36!=================================================================
37!
38!      This is a modified version that accepts spectrally varying input
39!      RW (2010)
40!
41!=================================================================
42 
43! Addition by RW (2010) to allow OSR to be saved in .nc format
44      use radinc_h, only : L_NSPECTV
45      use geometry_mod, only: cell_area
46      use mod_phys_lmdz_para, only : is_mpi_root, is_master, gather
47      use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo,
48     &                              nbp_lon, nbp_lat, grid_type,
49     &                              unstructured
50      use time_phylmdz_mod, only: ecritphy, iphysiq, day_step, day_ini
51      use callkeys_mod, only: iradia
52
53      implicit none
54
55      include "netcdf.inc"
56
57! Arguments on input:
58      integer ngrid
59      character (len=*) :: nom,titre,unite
60      integer dimpx
61      real px(ngrid,L_NSPECTV)
62
63! Local variables:
64
65!      real dx3(iip1,jjp1,llm) ! to store a 3D data set
66!      real dx2(iip1,jjp1)     ! to store a 2D (surface) data set
67!      real dx0
68
69      integer irythme
70      integer ierr
71      integer iq
72      integer i,j,l,zmax , ig0
73
74      integer,save :: zitau=0
75      character(len=20),save :: firstnom='1234567890'
76      real,save :: date
77!$OMP THREADPRIVATE(firstnom,zitau,date)
78
79! Ajouts
80      integer, save :: ntime=0
81!$OMP THREADPRIVATE(ntime)
82      integer :: idim,varid
83      integer :: nid
84      character (len =50):: fichnom
85      integer, dimension(4) :: id
86      integer, dimension(4) :: edges,corner
87
88      real area((nbp_lon+1),nbp_lat)
89! added by RDW for OSR output
90      real dx3(nbp_lon+1,nbp_lat,L_NSPECTV) ! to store the data set
91      real dx3_1d(1,L_NSPECTV) ! to store the data with 1D model
92
93#ifdef CPP_PARA
94! Added to work in parallel mode
95      real dx3_glop(klon_glo,L_NSPECTV)
96      real dx3_glo(nbp_lon,nbp_lat,L_NSPECTV) ! to store a global 3D data set
97      real areafi_glo(klon_glo) ! mesh area on global physics grid
98#else
99      real areafi_glo(ngrid) ! mesh area on global physics grid
100#endif
101      if (grid_type == unstructured) then
102        return
103      endif
104
105!***************************************************************
106!Sortie des variables au rythme voulu
107
108      irythme = ecritphy*iradia ! sortie au rythme de ecritphy
109!EM+JL if the spetra need to be output more frequently, need to define a ecritSpec...
110!     irythme = iphysiq  ! sortie a tous les pas physique
111
112!***************************************************************
113
114! Initialisation of 'firstnom' and create/open the "diagfi.nc" NetCDF file
115! ------------------------------------------------------------------------
116! (Au tout premier appel de la subroutine durant le run.)
117
118      fichnom="diagspecVI.nc"
119
120      if (firstnom.eq.'1234567890') then ! .true. for the very first call
121      !  to this subroutine; now set 'firstnom'
122         firstnom = nom
123         ! just to be sure, check that firstnom is large enough to hold nom
124         if (len_trim(firstnom).lt.len_trim(nom)) then
125           write(*,*) "writediagfi: Error !!!"
126           write(*,*) "   firstnom string not long enough!!"
127           write(*,*) "   increase its size to at least ",len_trim(nom)
128           stop
129         endif
130
131#ifdef CPP_PARA
132          ! Gather cell_area() mesh area on physics grid
133          call Gather(cell_area,areafi_glo)
134#else
135         areafi_glo(:)=cell_area(:)
136#endif
137         ! Create the NetCDF file
138         if (is_master) then
139         ierr = NF_CREATE(fichnom, NF_CLOBBER, nid)
140         ! Define the 'Time' dimension
141         ierr = nf_def_dim(nid,"Time",NF_UNLIMITED,idim)
142         ! Define the 'Time' variable
143#ifdef NC_DOUBLE
144         ierr = NF_DEF_VAR (nid, "Time", NF_DOUBLE, 1, idim,varid)
145#else
146         ierr = NF_DEF_VAR (nid, "Time", NF_FLOAT, 1, idim,varid)
147#endif
148         ! Add a long_name attribute
149         ierr = NF_PUT_ATT_TEXT (nid, varid, "long_name",
150     .          4,"Time")
151         ! Add a units attribute
152         ierr = NF_PUT_ATT_TEXT(nid, varid,'units',29,
153     .          "days since 0000-00-0 00:00:00")
154         ! Switch out of NetCDF Define mode
155         ierr = NF_ENDDEF(nid)
156
157         ! Build area()
158         IF (klon_glo>1) THEN
159          do i=1,nbp_lon+1 ! poles
160           ! divide at the poles by nbp_lon
161           area(i,1)=areafi_glo(1)/nbp_lon
162           area(i,nbp_lat)=areafi_glo(klon_glo)/nbp_lon
163          enddo
164          do j=2,nbp_lat-1
165           ig0= 1+(j-2)*nbp_lon
166           do i=1,nbp_lon
167              area(i,j)=areafi_glo(ig0+i)
168           enddo
169           ! handle redundant point in longitude
170           area(nbp_lon+1,j)=area(1,j)
171          enddo
172         ENDIF
173
174         ! write "header" of file (longitudes, latitudes, geopotential, ...)
175         IF (klon_glo>1) THEN ! general 3D case
176           call iniwrite_specVI(nid,day_ini,area,nbp_lon+1,nbp_lat)
177         ELSE
178           call iniwrite_specVI(nid,day_ini,areafi_glo(1),1,1)
179         ENDIF
180         endif ! of if (is_master)
181
182         zitau = -1 ! initialize zitau
183      else
184         if (is_master) then
185           ! Open the NetCDF file
186           ierr = NF_OPEN(fichnom,NF_WRITE,nid)
187         endif
188      endif ! if (firstnom.eq.'1234567890')
189
190! Increment time index 'zitau' if it is the "firstcall" (at given time level)
191! to writediagfi
192!------------------------------------------------------------------------
193      if (nom.eq.firstnom) then
194          zitau = zitau + iphysiq
195      end if
196
197!--------------------------------------------------------
198! Write the variables to output file if it's time to do so
199!--------------------------------------------------------
200
201      if ( MOD(zitau+1,irythme) .eq.0.) then
202
203! Compute/write/extend 'Time' coordinate (date given in days)
204! (done every "first call" (at given time level) to writediagfi)
205! Note: date is incremented as 1 step ahead of physics time
206!       (like the 'histoire' outputs)
207!--------------------------------------------------------
208
209        if (nom.eq.firstnom) then
210
211        ! We have identified a "first call" (at given date)
212           ntime=ntime+1 ! increment # of stored time steps
213           ! compute corresponding date (in days and fractions thereof)
214           date= float (zitau +1)/float (day_step)
215
216           if (is_master) then
217             ! Get NetCDF ID of 'Time' variable
218             ierr= NF_INQ_VARID(nid,"Time",varid)
219
220             ! Write (append) the new date to the 'Time' array
221#ifdef NC_DOUBLE
222             ierr= NF_PUT_VARA_DOUBLE(nid,varid,ntime,1,date)
223#else
224             ierr= NF_PUT_VARA_REAL(nid,varid,ntime,1,date)
225#endif
226             if (ierr.ne.NF_NOERR) then
227              write(*,*) "***** PUT_VAR matter in writediagspec_nc"
228              write(*,*) "***** with time"
229              write(*,*) 'ierr=', ierr   
230c             call abort
231             endif
232
233             write(6,*)'WRITEDIAGSPECVI: date= ', date
234           endif ! of if (is_master)
235        end if ! of if (nom.eq.firstnom)
236
237
238 
239!Case of a 3D variable
240!---------------------
241        if (dimpx.eq.3) then
242
243!         A. Recast (copy) variable from physics grid to dynamics grid
244#ifdef CPP_PARA
245  ! gather field on a "global" (without redundant longitude) array
246          call Gather(px,dx3_glop)
247!$OMP MASTER
248          if (is_mpi_root) then
249            call Grid1Dto2D_glo(dx3_glop,dx3_glo)
250            ! copy dx3_glo() to dx3(:) and add redundant longitude
251            dx3(1:nbp_lon,:,:)=dx3_glo(1:nbp_lon,:,:)
252            dx3(nbp_lon+1,:,:)=dx3(1,:,:)
253          endif
254!$OMP END MASTER
255!$OMP BARRIER
256#else
257          IF (klon_glo>1) THEN ! General case
258           DO l=1,L_NSPECTV
259             DO i=1,nbp_lon+1
260                dx3(i,1,l)=px(1,l)
261                dx3(i,nbp_lat,l)=px(ngrid,l)
262             ENDDO
263             DO j=2,nbp_lat-1
264                ig0= 1+(j-2)*nbp_lon
265                DO i=1,nbp_lon
266                   dx3(i,j,l)=px(ig0+i,l)
267                ENDDO
268                dx3(nbp_lon+1,j,l)=dx3(1,j,l)
269             ENDDO
270           ENDDO
271          ELSE ! 1D model case
272           dx3_1d(1,1:L_NSPECTV)=px(1,1:L_NSPECTV)
273          ENDIF
274#endif
275
276!         B. Write (append) the variable to the NetCDF file
277          if (is_master) then
278
279! name of the variable
280           ierr= NF_INQ_VARID(nid,nom,varid)
281           if (ierr /= NF_NOERR) then
282! corresponding dimensions
283              ierr= NF_INQ_DIMID(nid,"longitude",id(1))
284              ierr= NF_INQ_DIMID(nid,"latitude",id(2))
285              ierr= NF_INQ_DIMID(nid,"VI_Wavenumber",id(3))
286              ierr= NF_INQ_DIMID(nid,"Time",id(4))
287
288! Create the variable if it doesn't exist yet
289
290              write (*,*) "=========================="
291              write (*,*) "DIAGSPEC: creating variable ",nom
292              call def_var(nid,nom,titre,unite,4,id,varid,ierr)
293
294           endif
295
296           corner(1)=1
297           corner(2)=1
298           corner(3)=1
299           corner(4)=ntime
300
301           IF (klon_glo==1) THEN
302             edges(1)=1
303           ELSE
304             edges(1)=nbp_lon+1
305           ENDIF
306           edges(2)=nbp_lat
307           edges(3)=L_NSPECTV
308           edges(4)=1
309#ifdef NC_DOUBLE
310           IF (klon_glo>1) THEN ! General case
311             ierr= NF_PUT_VARA_DOUBLE(nid,varid,corner,edges,dx3)
312           ELSE
313             ierr= NF_PUT_VARA_DOUBLE(nid,varid,corner,edges,dx3_1d)
314           ENDIF
315#else
316           IF (klon_glo>1) THEN ! General case
317             ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx3)
318           ELSE
319             ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx3_1d)
320           ENDIF
321#endif
322
323           if (ierr.ne.NF_NOERR) then
324              write(*,*) "***** PUT_VAR problem in writediagspec"
325              write(*,*) "***** with ",nom
326              write(*,*) 'ierr=', ierr
327             call abort
328           endif
329
330          endif ! of if (is_master)
331
332        endif ! of if (dimpx.eq.3)
333
334      endif ! of if ( MOD(zitau+1,irythme) .eq.0.)
335
336      ! Close the NetCDF file
337      if (is_master) then
338        ierr= NF_CLOSE(nid)
339      endif
340
341      end
Note: See TracBrowser for help on using the repository browser.