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

Last change on this file since 3562 was 2732, checked in by aslmd, 3 years ago

took into account unstructured grid in writediagspec*.F

File size: 11.3 KB
RevLine 
[965]1      subroutine writediagspecVI(ngrid,nom,titre,unite,dimpx,px)
[305]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)
[965]34!      dimpx : dimension de px : 0, 2, ou 3 dimensions
[305]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
[1543]45      use geometry_mod, only: cell_area
[965]46      use mod_phys_lmdz_para, only : is_mpi_root, is_master, gather
[1529]47      use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo,
[2732]48     &                              nbp_lon, nbp_lat, grid_type,
49     &                              unstructured
[1525]50      use time_phylmdz_mod, only: ecritphy, iphysiq, day_step, day_ini
[1397]51      use callkeys_mod, only: iradia
[305]52
53      implicit none
54
[1529]55      include "netcdf.inc"
[305]56
57! Arguments on input:
58      integer ngrid
59      character (len=*) :: nom,titre,unite
[965]60      integer dimpx
[305]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
[1529]74      integer,save :: zitau=0
75      character(len=20),save :: firstnom='1234567890'
76      real,save :: date
[1315]77!$OMP THREADPRIVATE(firstnom,zitau,date)
[305]78
79! Ajouts
80      integer, save :: ntime=0
[1315]81!$OMP THREADPRIVATE(ntime)
[305]82      integer :: idim,varid
83      integer :: nid
84      character (len =50):: fichnom
85      integer, dimension(4) :: id
86      integer, dimension(4) :: edges,corner
87
[1529]88      real area((nbp_lon+1),nbp_lat)
[305]89! added by RDW for OSR output
[1531]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
[305]92
[965]93#ifdef CPP_PARA
94! Added to work in parallel mode
95      real dx3_glop(klon_glo,L_NSPECTV)
[1529]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
[965]98#else
[1529]99      real areafi_glo(ngrid) ! mesh area on global physics grid
[965]100#endif
[2732]101      if (grid_type == unstructured) then
102        return
103      endif
[305]104
105!***************************************************************
106!Sortie des variables au rythme voulu
107
[1216]108      irythme = ecritphy*iradia ! sortie au rythme de ecritphy
[533]109!EM+JL if the spetra need to be output more frequently, need to define a ecritSpec...
[305]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
[1529]131#ifdef CPP_PARA
[1542]132          ! Gather cell_area() mesh area on physics grid
133          call Gather(cell_area,areafi_glo)
[1529]134#else
[1542]135         areafi_glo(:)=cell_area(:)
[1529]136#endif
[305]137         ! Create the NetCDF file
[965]138         if (is_master) then
[305]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
[1529]157         ! Build area()
[1531]158         IF (klon_glo>1) THEN
159          do i=1,nbp_lon+1 ! poles
[1529]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
[1531]163          enddo
164          do j=2,nbp_lat-1
[1529]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)
[1531]171          enddo
172         ENDIF
[1529]173
[305]174         ! write "header" of file (longitudes, latitudes, geopotential, ...)
[1531]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
[965]180         endif ! of if (is_master)
[305]181
182         zitau = -1 ! initialize zitau
183      else
[965]184         if (is_master) then
185           ! Open the NetCDF file
186           ierr = NF_OPEN(fichnom,NF_WRITE,nid)
187         endif
[305]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
[965]216           if (is_master) then
217             ! Get NetCDF ID of 'Time' variable
218             ierr= NF_INQ_VARID(nid,"Time",varid)
[305]219
[965]220             ! Write (append) the new date to the 'Time' array
[305]221#ifdef NC_DOUBLE
[965]222             ierr= NF_PUT_VARA_DOUBLE(nid,varid,ntime,1,date)
[305]223#else
[965]224             ierr= NF_PUT_VARA_REAL(nid,varid,ntime,1,date)
[305]225#endif
[965]226             if (ierr.ne.NF_NOERR) then
[305]227              write(*,*) "***** PUT_VAR matter in writediagspec_nc"
228              write(*,*) "***** with time"
229              write(*,*) 'ierr=', ierr   
230c             call abort
[965]231             endif
[305]232
[1529]233             write(6,*)'WRITEDIAGSPECVI: date= ', date
[965]234           endif ! of if (is_master)
[305]235        end if ! of if (nom.eq.firstnom)
236
237
238 
239!Case of a 3D variable
240!---------------------
[965]241        if (dimpx.eq.3) then
[305]242
[965]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
[1529]251            dx3(1:nbp_lon,:,:)=dx3_glo(1:nbp_lon,:,:)
252            dx3(nbp_lon+1,:,:)=dx3(1,:,:)
[965]253          endif
254!$OMP END MASTER
255!$OMP BARRIER
256#else
[1531]257          IF (klon_glo>1) THEN ! General case
[305]258           DO l=1,L_NSPECTV
[1529]259             DO i=1,nbp_lon+1
[305]260                dx3(i,1,l)=px(1,l)
[1529]261                dx3(i,nbp_lat,l)=px(ngrid,l)
[305]262             ENDDO
[1529]263             DO j=2,nbp_lat-1
264                ig0= 1+(j-2)*nbp_lon
265                DO i=1,nbp_lon
[305]266                   dx3(i,j,l)=px(ig0+i,l)
267                ENDDO
[1529]268                dx3(nbp_lon+1,j,l)=dx3(1,j,l)
[305]269             ENDDO
270           ENDDO
[1531]271          ELSE ! 1D model case
272           dx3_1d(1,1:L_NSPECTV)=px(1,1:L_NSPECTV)
273          ENDIF
[965]274#endif
[305]275
[965]276!         B. Write (append) the variable to the NetCDF file
277          if (is_master) then
278
[305]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))
[1531]285              ierr= NF_INQ_DIMID(nid,"VI_Wavenumber",id(3))
[305]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
[1531]301           IF (klon_glo==1) THEN
302             edges(1)=1
303           ELSE
304             edges(1)=nbp_lon+1
305           ENDIF
[1529]306           edges(2)=nbp_lat
[305]307           edges(3)=L_NSPECTV
308           edges(4)=1
309#ifdef NC_DOUBLE
[1531]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
[305]315#else
[1531]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
[305]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
[965]330          endif ! of if (is_master)
[305]331
[965]332        endif ! of if (dimpx.eq.3)
333
[305]334      endif ! of if ( MOD(zitau+1,irythme) .eq.0.)
335
[965]336      ! Close the NetCDF file
337      if (is_master) then
338        ierr= NF_CLOSE(nid)
339      endif
[305]340
341      end
Note: See TracBrowser for help on using the repository browser.