source: IOIPSL/trunk/src/histcom.f90 @ 358

Last change on this file since 358 was 358, checked in by bellier, 14 years ago

A little step towards the CF Metadata Convention

  • Property svn:keywords set to Id
File size: 83.8 KB
Line 
1MODULE histcom
2!-
3!$Id$
4!-
5  USE netcdf
6!-
7  USE stringop,  ONLY : nocomma,cmpblank,findpos,find_str,strlowercase
8  USE mathelp,   ONLY : mathop,moycum,trans_buff,buildop
9  USE fliocom,   ONLY : flio_dom_file,flio_dom_att
10  USE calendar
11  USE errioipsl, ONLY : ipslerr
12!-
13  IMPLICIT NONE
14!-
15  PRIVATE
16  PUBLIC :: histbeg, histdef, histhori, histvert, histend, &
17 &          histwrite, histclo, histsync, ioconf_modname
18!---------------------------------------------------------------------
19!- Some confusing vocabulary in this code !
20!- =========================================
21!-
22!- A REGULAR grid is a grid which is i,j indexes
23!- and thus it is stored in a 2D matrix.
24!- This is opposed to a IRREGULAR grid which is only in a vector
25!- and where we do not know which neighbors we have.
26!- As a consequence we need the bounds for each grid-cell.
27!-
28!- A RECTILINEAR grid is a special case of a regular grid
29!- in which all longitudes for i constant are equal
30!- and all latitudes for j constant.
31!- In other words we do not need the full 2D matrix
32!- to describe the grid, just two vectors.
33!---------------------------------------------------------------------
34  INTERFACE histwrite
35!---------------------------------------------------------------------
36!- The "histwrite" routines will give the data to the I/O system.
37!- It will trigger the operations to be performed,
38!- and the writting to the file if needed
39!-
40!- We test for the work to be done at this time here so that at a
41!- later stage we can call different operation and write subroutine
42!- for the REAL and INTEGER interfaces
43!-
44!- INPUT
45!- pfileid  : The ID of the file on which this variable is to be,
46!-            written. The variable should have been defined in
47!-            this file before.
48!- pvarname : The short name of the variable
49!- pitau    : Current timestep
50!- pdata    : The variable, I mean the real data !
51!- nbindex  : The number of indexes provided. If it is equal to
52!-            the size of the full field as provided in histdef
53!-            then nothing is done.
54!- nindex   : The indices used to expand the variable (pdata)
55!-            onto the full field.
56!---------------------------------------------------------------------
57!- histwrite - we have to prepare different type of fields :
58!-             real and integer, 1,2 or 3D
59    MODULE PROCEDURE histwrite_r1d,histwrite_r2d,histwrite_r3d
60  END INTERFACE
61!-
62  INTERFACE histbeg
63    MODULE PROCEDURE histbeg_totreg,histbeg_regular,histbeg_irregular
64  END INTERFACE
65!-
66  INTERFACE histhori
67    MODULE PROCEDURE histhori_regular,histhori_irregular
68  END INTERFACE
69!-
70! Fixed parameter
71!-
72  INTEGER,PARAMETER :: nb_files_max=20, nb_var_max=400, &
73 &                     nb_hax_max=5, nb_zax_max=10, nbopp_max=10
74  REAL,PARAMETER :: missing_val=nf90_fill_real
75!-
76  INTEGER :: bufftmp_max(nb_files_max) = 1
77!-
78! Time variables
79!-
80  INTEGER,SAVE :: itau0(nb_files_max)=0
81  REAL,DIMENSION(nb_files_max),SAVE ::date0,deltat
82!-
83! Counter of elements
84!-
85  INTEGER,SAVE :: nb_files=0
86  INTEGER,DIMENSION(nb_files_max),SAVE :: nb_var=0,nb_tax=0
87!-
88! DOMAIN IDs for files
89!-
90  INTEGER,DIMENSION(nb_files_max),SAVE :: dom_id_svg=-1
91!-
92! NETCDF IDs for files and axes
93!-
94  INTEGER,DIMENSION(nb_files_max),SAVE :: ncdf_ids,xid,yid,tid
95!-
96! General definitions in the NETCDF file
97!-
98  INTEGER,DIMENSION(nb_files_max,2),SAVE :: &
99 &   full_size=0,slab_ori,slab_sz
100!-
101! The horizontal axes
102!-
103  INTEGER,SAVE :: nb_hax(nb_files_max)=0
104  CHARACTER(LEN=25),SAVE :: hax_name(nb_files_max,nb_hax_max,2)
105!-
106! The vertical axes
107!-
108  INTEGER,SAVE :: nb_zax(nb_files_max)=0
109  INTEGER,DIMENSION(nb_files_max,nb_zax_max),SAVE :: &
110  &  zax_size,zax_ids
111  CHARACTER(LEN=20),SAVE :: zax_name(nb_files_max,nb_zax_max)
112!-
113! Informations on each variable
114!-
115  INTEGER,DIMENSION(nb_files_max,nb_var_max),SAVE :: &
116 &  nbopp
117  CHARACTER(LEN=20),DIMENSION(nb_files_max,nb_var_max),SAVE :: &
118 &  name,unit_name
119  CHARACTER(LEN=80),DIMENSION(nb_files_max,nb_var_max),SAVE :: &
120 &  title,fullop
121  CHARACTER(LEN=7),SAVE :: topp(nb_files_max,nb_var_max)
122  CHARACTER(LEN=7),SAVE :: sopps(nb_files_max,nb_var_max,nbopp_max)
123  REAL,SAVE :: scal(nb_files_max,nb_var_max,nbopp_max)
124!- Sizes of the associated grid and zommed area
125  INTEGER,DIMENSION(nb_files_max,nb_var_max,3),SAVE :: &
126 &   scsize,zorig,zsize
127!- Sizes for the data as it goes through the various math operations
128  INTEGER,SAVE :: datasz_in(nb_files_max,nb_var_max,3) = -1
129  INTEGER,SAVE :: datasz_max(nb_files_max,nb_var_max) = -1
130!-
131  INTEGER,DIMENSION(nb_files_max,nb_var_max),SAVE :: &
132 &  var_haxid,var_zaxid,var_axid,ncvar_ids
133!-
134  REAL,DIMENSION(nb_files_max,nb_var_max,2),SAVE :: hist_minmax
135  LOGICAL,DIMENSION(nb_files_max,nb_var_max),SAVE :: &
136 &  hist_calc_rng=.FALSE.,hist_wrt_rng=.FALSE.
137!-
138  REAL,DIMENSION(nb_files_max,nb_var_max),SAVE :: &
139 &  freq_opp,freq_wrt
140  INTEGER,DIMENSION(nb_files_max,nb_var_max),SAVE :: &
141 &  last_opp,last_wrt,last_opp_chk,last_wrt_chk,nb_opp,nb_wrt,point
142!-
143! Book keeping for the buffers
144!-
145  INTEGER,SAVE :: buff_pos=0
146  REAL,ALLOCATABLE,SAVE :: buffer(:)
147  LOGICAL,SAVE :: &
148 &  zoom(nb_files_max)=.FALSE., regular(nb_files_max)=.TRUE.
149!-
150! Book keeping of the axes
151!-
152  INTEGER,DIMENSION(nb_files_max,nb_var_max),SAVE :: &
153 &  tdimid,tax_last
154  CHARACTER(LEN=40),DIMENSION(nb_files_max,nb_var_max),SAVE :: &
155 &  tax_name
156!-
157! A list of functions which require special action
158! (Needs to be updated when functions are added
159!  but they are well located here)
160!-
161  CHARACTER(LEN=120),SAVE :: &
162 &  indchfun = 'scatter, fill, gather, coll', &
163 &  fuchnbout = 'scatter, fill'
164!- Some configurable variables with locks
165  CHARACTER(LEN=80),SAVE :: model_name='An IPSL model'
166  LOGICAL,SAVE :: lock_modname=.FALSE.
167!-
168!===
169CONTAINS
170!===
171!-
172SUBROUTINE histbeg_totreg                     &
173 & (pfilename, pim, plon, pjm, plat,          &
174 &  par_orix, par_szx, par_oriy, par_szy,     &
175 &  pitau0, pdate0, pdeltat, phoriid, pfileid, domain_id)
176!---------------------------------------------------------------------
177!- This is just an interface for histbeg_regular in case when
178!- the user provides plon and plat as vectors.
179!- Obviously this can only be used for very regular grids.
180!-
181!- INPUT
182!-
183!- pfilename : Name of the netcdf file to be created
184!- pim       : Size of arrays in longitude direction
185!- plon      : Coordinates of points in longitude
186!- pjm       : Size of arrays in latitude direction
187!- plat      : Coordinates of points in latitude
188!-
189!- The next 4 arguments allow to define a horizontal zoom
190!- for this file. It is assumed that all variables to come
191!- have the same index space. This can not be assumed for
192!- the z axis and thus we define the zoom in histdef.
193!-
194!- par_orix  : Origin of the slab of data within the X axis (pim)
195!- par_szx   : Size of the slab of data in X
196!- par_oriy  : Origin of the slab of data within the Y axis (pjm)
197!- par_szy   : Size of the slab of data in Y
198!-
199!- pitau0    : time step at which the history tape starts
200!- pdate0    : The Julian date at which the itau was equal to 0
201!- pdeltat   : Time step in seconds. Time step of the counter itau
202!-             used in histwrt for instance
203!-
204!- OUTPUT
205!-
206!- phoriid   : ID of the horizontal grid
207!- pfileid   : ID of the netcdf file
208!-
209!- Optional INPUT arguments
210!-
211!- domain_id       : Domain identifier
212!-
213!- TO DO
214!-
215!- This package should be written in f90
216!- and use the following features :
217!-    - structures for the meta-data of the files and variables
218!-    - memory allocation as needed
219!-    - Pointers
220!-
221!- VERSION
222!-
223!---------------------------------------------------------------------
224  IMPLICIT NONE
225!-
226  CHARACTER(LEN=*) :: pfilename
227  INTEGER,INTENT(IN) :: pim,pjm
228  REAL,DIMENSION(pim),INTENT(IN) :: plon
229  REAL,DIMENSION(pjm),INTENT(IN) :: plat
230  INTEGER,INTENT(IN):: par_orix,par_szx,par_oriy,par_szy
231  INTEGER,INTENT(IN) :: pitau0
232  REAL,INTENT(IN) :: pdate0,pdeltat
233  INTEGER,INTENT(OUT) :: pfileid,phoriid
234  INTEGER,INTENT(IN),OPTIONAL :: domain_id
235!-
236  LOGICAL :: check = .FALSE.
237!-
238  REAL,ALLOCATABLE,DIMENSION(:,:) :: lon_tmp,lat_tmp
239!---------------------------------------------------------------------
240  IF (check) WRITE(*,*) "histbeg_totreg"
241!-
242  ALLOCATE (lon_tmp(pim,pjm),lat_tmp(pim,pjm))
243!-
244  lon_tmp(:,:) = SPREAD(plon(:),2,pjm)
245  lat_tmp(:,:) = SPREAD(plat(:),1,pim)
246!-
247  CALL histbeg_regular &
248 &  (pfilename,pim,lon_tmp,pjm,lat_tmp, &
249 &   par_orix,par_szx,par_oriy,par_szy, &
250 &   pitau0,pdate0,pdeltat,phoriid,pfileid, &
251 &   .TRUE.,domain_id)
252!-
253  DEALLOCATE (lon_tmp, lat_tmp)
254!----------------------------
255END SUBROUTINE histbeg_totreg
256!===
257SUBROUTINE histbeg_regular &
258 & (pfilename, pim, plon, pjm, plat,         &
259 &  par_orix, par_szx, par_oriy, par_szy,    &
260 &  pitau0, pdate0, pdeltat, phoriid, pfileid, &
261 &  opt_rectilinear, domain_id)
262!---------------------------------------------------------------------
263!- This subroutine initializes a netcdf file and returns the ID.
264!- It will set up the geographical space on which the data will be
265!- stored and offers the possibility of seting a zoom.
266!- It also gets the global parameters into the I/O subsystem.
267!-
268!- INPUT
269!-
270!- pfilename : Name of the netcdf file to be created
271!- pim       : Size of arrays in longitude direction
272!- plon      : Coordinates of points in longitude
273!- pjm       : Size of arrays in latitude direction
274!- plat      : Coordinates of points in latitude
275!-
276!- The next 4 arguments allow to define a horizontal zoom
277!- for this file. It is assumed that all variables to come
278!- have the same index space. This can not be assumed for
279!- the z axis and thus we define the zoom in histdef.
280!-
281!- par_orix  : Origin of the slab of data within the X axis (pim)
282!- par_szx   : Size of the slab of data in X
283!- par_oriy  : Origin of the slab of data within the Y axis (pjm)
284!- par_szy   : Size of the slab of data in Y
285!-
286!- pitau0    : time step at which the history tape starts
287!- pdate0    : The Julian date at which the itau was equal to 0
288!- pdeltat   : Time step in seconds. Time step of the counter itau
289!-             used in histwrt for instance
290!-
291!- OUTPUT
292!-
293!- phoriid   : ID of the horizontal grid
294!- pfileid   : ID of the netcdf file
295!-
296!- Optional INPUT arguments
297!-
298!- opt_rectilinear : If true we know the grid is rectilinear
299!- domain_id       : Domain identifier
300!-
301!- TO DO
302!-
303!- This package should be written in F90 and use the following
304!- feature :
305!-   - structures for the meta-data of the files and variables
306!-   - memory allocation as needed
307!-   - Pointers
308!-
309!- VERSION
310!-
311!---------------------------------------------------------------------
312  IMPLICIT NONE
313!-
314  CHARACTER(LEN=*) :: pfilename
315  INTEGER,INTENT(IN) :: pim,pjm
316  REAL,DIMENSION(pim,pjm),INTENT(IN) :: plon,plat
317  INTEGER,INTENT(IN):: par_orix,par_szx,par_oriy,par_szy
318  INTEGER,INTENT(IN) :: pitau0
319  REAL,INTENT(IN) :: pdate0,pdeltat
320  INTEGER,INTENT(OUT) :: pfileid,phoriid
321  LOGICAL,INTENT(IN),OPTIONAL :: opt_rectilinear
322  INTEGER,INTENT(IN),OPTIONAL :: domain_id
323!-
324  INTEGER :: ncid,iret
325  CHARACTER(LEN=120) :: file
326  CHARACTER(LEN=30) :: timenow
327  LOGICAL :: rectilinear
328!-
329  LOGICAL :: check = .FALSE.
330!---------------------------------------------------------------------
331  nb_files = nb_files+1
332  pfileid  = nb_files
333!-
334! 1.0 Transfering into the common for future use
335!-
336  IF (check) WRITE(*,*) "histbeg_regular 1.0"
337!-
338  itau0(pfileid) = pitau0
339  date0(pfileid) = pdate0
340  deltat(pfileid) = pdeltat
341!-
342  IF (PRESENT(opt_rectilinear)) THEN
343    rectilinear = opt_rectilinear
344  ELSE
345    rectilinear = .FALSE.
346  ENDIF
347!-
348! 2.0 Initializes all variables for this file
349!-
350  IF (check) WRITE(*,*) "histbeg_regular 2.0"
351!-
352  IF (nb_files > nb_files_max) THEN
353    CALL ipslerr (3,"histbeg", &
354   &  'Table of files too small. You should increase nb_files_max', &
355   &  'in histcom.f90 in order to accomodate all these files', ' ')
356  ENDIF
357!-
358  nb_var(pfileid) = 0
359  nb_tax(pfileid) = 0
360  nb_hax(pfileid) = 0
361  nb_zax(pfileid) = 0
362!-
363  slab_ori(pfileid,1:2) = (/ par_orix, par_oriy /)
364  slab_sz(pfileid,1:2)  = (/ par_szx,  par_szy  /)
365!-
366! 3.0 Opening netcdf file and defining dimensions
367!-
368  IF (check) WRITE(*,*) "histbeg_regular 3.0"
369!-
370! Add DOMAIN number and ".nc" suffix in file name if needed
371!-
372  file = pfilename
373  CALL flio_dom_file (file,domain_id)
374!-
375  iret = NF90_CREATE (file,NF90_CLOBBER,ncid)
376!-
377  IF (rectilinear) THEN
378    iret = NF90_DEF_DIM (ncid,'lon',par_szx,xid(nb_files))
379    iret = NF90_DEF_DIM (ncid,'lat',par_szy,yid(nb_files))
380  ELSE
381    iret = NF90_DEF_DIM (ncid,'x',par_szx,xid(nb_files))
382    iret = NF90_DEF_DIM (ncid,'y',par_szy,yid(nb_files))
383  ENDIF
384!-
385! 4.0 Declaring the geographical coordinates and other attributes
386!-
387  IF (check) WRITE(*,*) "histbeg_regular 4.0"
388!-
389! 4.3 Global attributes
390!-
391  iret = NF90_PUT_ATT (ncid,NF90_GLOBAL,'Conventions','CF-1.3')
392  iret = NF90_PUT_ATT (ncid,NF90_GLOBAL,'file_name',TRIM(file))
393  iret = NF90_PUT_ATT (ncid,NF90_GLOBAL,'production',TRIM(model_name))
394  lock_modname = .TRUE.
395  CALL ioget_timestamp (timenow)
396  iret = NF90_PUT_ATT (ncid,NF90_GLOBAL,'TimeStamp',TRIM(timenow))
397!-
398! 5.0 Saving some important information on this file in the common
399!-
400  IF (check) WRITE(*,*) "histbeg_regular 5.0"
401!-
402  IF (PRESENT(domain_id)) THEN
403    dom_id_svg(pfileid) = domain_id
404  ENDIF
405  ncdf_ids(pfileid) = ncid
406  full_size(pfileid,1:2) = (/ pim, pjm /)
407!-
408! 6.0 storing the geographical coordinates
409!-
410  zoom(pfileid) = (pim /= par_szx).OR.(pjm /= par_szy)
411  regular(pfileid)=.TRUE.
412!-
413  CALL histhori_regular (pfileid, pim, plon, pjm, plat, &
414 &  ' ' , 'Default grid', phoriid, rectilinear)
415!-----------------------------
416END SUBROUTINE histbeg_regular
417!===
418SUBROUTINE histbeg_irregular &
419 &  (pfilename, pim, plon, plon_bounds, plat, plat_bounds,   &
420 &   pitau0, pdate0, pdeltat, phoriid, pfileid, domain_id)
421!---------------------------------------------------------------------
422!- This subroutine initializes a netcdf file and returns the ID.
423!- This version is for totaly irregular grids. In this case all
424!- all the data comes in as vectors and for the grid we have
425!- the coordinates of the 4 corners.
426!- It also gets the global parameters into the I/O subsystem.
427!-
428!- INPUT
429!-
430!- pfilename   : Name of the netcdf file to be created
431!- pim         : Size of arrays in longitude direction
432!- plon        : Coordinates of points in longitude
433!- plon_bounds : The 2 corners of the grid in longitude
434!- plat        : Coordinates of points in latitude
435!- plat_bounds : The 2 corners of the grid in latitude
436!-
437!- pitau0    : time step at which the history tape starts
438!- pdate0    : The Julian date at which the itau was equal to 0
439!- pdeltat   : Time step in seconds. Time step of the counter itau
440!-             used in histwrt for instance
441!-
442!- OUTPUT
443!-
444!- phoriid   : ID of the horizontal grid
445!- pfileid   : ID of the netcdf file
446!-
447!- Optional INPUT arguments
448!-
449!- domain_id       : Domain identifier
450!-
451!- TO DO
452!-
453!- This package should be written in F90 and use the following
454!- feature :
455!-   - structures for the meta-data of the files and variables
456!-   - memory allocation as needed
457!-   - Pointers
458!-
459!- VERSION
460!-
461!---------------------------------------------------------------------
462  IMPLICIT NONE
463!-
464  CHARACTER(LEN=*) :: pfilename
465  INTEGER,INTENT(IN) :: pim
466  REAL,DIMENSION(pim),INTENT(IN) :: plon,plat
467  REAL,DIMENSION(:,:),INTENT(IN) :: plon_bounds,plat_bounds
468  INTEGER,INTENT(IN) :: pitau0
469  REAL,INTENT(IN) :: pdate0,pdeltat
470  INTEGER,INTENT(OUT) :: pfileid,phoriid
471  INTEGER,INTENT(IN),OPTIONAL :: domain_id
472!-
473  INTEGER :: ncid,iret
474  CHARACTER(LEN=120) :: file
475  CHARACTER(LEN=30) :: timenow
476!-
477  LOGICAL :: check = .FALSE.
478!---------------------------------------------------------------------
479  nb_files = nb_files+1
480  pfileid  = nb_files
481!-
482! 1.0 Transfering into the common for future use
483!-
484  IF (check) WRITE(*,*) "histbeg_irregular 1.0"
485!-
486  itau0(pfileid) = pitau0
487  date0(pfileid) = pdate0
488  deltat(pfileid) = pdeltat
489!-
490! 2.0 Initializes all variables for this file
491!-
492  IF (check) WRITE(*,*) "histbeg_irregular 2.0"
493!-
494  IF (nb_files > nb_files_max) THEN
495    CALL ipslerr (3,"histbeg", &
496   &  'Table of files too small. You should increase nb_files_max', &
497   &  'in histcom.f90 in order to accomodate all these files',' ')
498  ENDIF
499!-
500  nb_var(pfileid) = 0
501  nb_tax(pfileid) = 0
502  nb_hax(pfileid) = 0
503  nb_zax(pfileid) = 0
504!-
505  slab_ori(pfileid,1:2) = (/ 1, 1 /)
506  slab_sz(pfileid,1:2)  = (/ pim, 1 /)
507!-
508! 3.0 Opening netcdf file and defining dimensions
509!-
510  IF (check) WRITE(*,*) "histbeg_irregular 3.0"
511!-
512! Add DOMAIN number and ".nc" suffix in file name if needed
513!-
514  file  = pfilename
515  CALL flio_dom_file (file,domain_id)
516!-
517  iret = NF90_CREATE (file,NF90_CLOBBER,ncid)
518!-
519  iret = NF90_DEF_DIM (ncid,'x',pim,xid(nb_files))
520  yid(nb_files) = 0
521!-
522!- 4.0 Declaring the geographical coordinates and other attributes
523!-
524   IF (check) WRITE(*,*) "histbeg_irregular 4.0"
525!-
526! 4.3 Global attributes
527!-
528  iret = NF90_PUT_ATT (ncid,NF90_GLOBAL,'Conventions','CF-1.3')
529  iret = NF90_PUT_ATT (ncid,NF90_GLOBAL,'file_name',TRIM(file))
530  iret = NF90_PUT_ATT (ncid,NF90_GLOBAL,'production',TRIM(model_name))
531  lock_modname = .TRUE.
532  CALL ioget_timestamp (timenow)
533  iret = NF90_PUT_ATT (ncid,NF90_GLOBAL,'TimeStamp',TRIM(timenow))
534!-
535! 5.0 Saving some important information on this file in the common
536!-
537  IF (check) WRITE(*,*) "histbeg_irregular 5.0"
538!-
539  IF (PRESENT(domain_id)) THEN
540    dom_id_svg(pfileid) = domain_id
541  ENDIF
542  ncdf_ids(pfileid) = ncid
543  full_size(pfileid,1:2) = (/ pim, 1 /)
544!-
545! 6.0 storing the geographical coordinates
546!-
547  zoom(pfileid)=.FALSE.
548  regular(pfileid)=.FALSE.
549!-
550  CALL histhori_irregular &
551 &  (pfileid, pim, plon, plon_bounds, plat, plat_bounds, &
552 &   ' ' , 'Default grid', phoriid)
553!-------------------------------
554END SUBROUTINE histbeg_irregular
555!===
556SUBROUTINE histhori_regular &
557 &  (pfileid,pim,plon,pjm,plat,phname,phtitle,phid,opt_rectilinear)
558!---------------------------------------------------------------------
559!- This subroutine is made to declare a new horizontale grid.
560!- It has to have the same number of points as
561!- the original and thus in this routine we will only
562!- add two variable (longitude and latitude).
563!- Any variable in the file can thus point to this pair
564!- through an attribute. This routine is very usefull
565!- to allow staggered grids.
566!-
567!- INPUT
568!-
569!- pfileid : The id of the file to which the grid should be added
570!- pim     : Size in the longitude direction
571!- plon    : The longitudes
572!- pjm     : Size in the latitude direction
573!- plat    : The latitudes
574!- phname  : The name of grid
575!- phtitle : The title of the grid
576!-
577!- OUTPUT
578!-
579!- phid    : Id of the created grid
580!-
581!- OPTIONAL
582!-
583!- opt_rectilinear : If true we know the grid is rectilinear.
584!-
585!---------------------------------------------------------------------
586  IMPLICIT NONE
587!-
588  INTEGER,INTENT(IN) :: pfileid, pim, pjm
589  REAL,INTENT(IN),DIMENSION(pim,pjm) :: plon, plat
590  CHARACTER(LEN=*),INTENT(IN) :: phname, phtitle
591  INTEGER,INTENT(OUT) :: phid
592  LOGICAL,INTENT(IN),OPTIONAL :: opt_rectilinear
593!-
594  CHARACTER(LEN=25) :: lon_name, lat_name
595  CHARACTER(LEN=80) :: tmp_title, tmp_name
596  INTEGER :: ndim
597  INTEGER,DIMENSION(2) :: dims
598  INTEGER :: nlonid, nlatid
599  INTEGER :: orix, oriy, par_szx, par_szy
600  INTEGER :: iret, ncid
601  LOGICAL :: rectilinear
602!-
603  LOGICAL :: check = .FALSE.
604!---------------------------------------------------------------------
605!-
606! 1.0 Check that all fits in the buffers
607!-
608  IF (    (pim /= full_size(pfileid,1)) &
609    & .OR.(pjm /= full_size(pfileid,2)) ) THEN
610    CALL ipslerr (3,"histhori", &
611   &  'The new horizontal grid does not have the same size', &
612   &  'as the one provided to histbeg. This is not yet ', &
613   &  'possible in the hist package.')
614  ENDIF
615!-
616  IF (PRESENT(opt_rectilinear)) THEN
617    rectilinear = opt_rectilinear
618  ELSE
619    rectilinear = .FALSE.
620  ENDIF
621!-
622! 1.1 Create all the variables needed
623!-
624  IF (check) WRITE(*,*) "histhori_regular 1.0"
625!-
626  ncid = ncdf_ids(pfileid)
627!-
628  ndim = 2
629  dims(1:2) = (/ xid(pfileid), yid(pfileid) /)
630!-
631  tmp_name = phname
632  IF (rectilinear) THEN
633    IF (nb_hax(pfileid) == 0) THEN
634      lon_name = 'lon'
635      lat_name = 'lat'
636    ELSE
637      lon_name = 'lon_'//TRIM(tmp_name)
638      lat_name = 'lat_'//TRIM(tmp_name)
639    ENDIF
640  ELSE
641    IF (nb_hax(pfileid) == 0) THEN
642      lon_name = 'nav_lon'
643      lat_name = 'nav_lat'
644    ELSE
645      lon_name = 'nav_lon_'//TRIM(tmp_name)
646      lat_name = 'nav_lat_'//TRIM(tmp_name)
647    ENDIF
648  ENDIF
649!-
650! 1.2 Save the informations
651!-
652  phid =  nb_hax(pfileid)+1
653  nb_hax(pfileid) = phid
654!-
655  hax_name(pfileid,phid,1:2) = (/ lon_name, lat_name /)
656  tmp_title = phtitle
657!-
658! 2.0 Longitude
659!-
660  IF (check) WRITE(*,*) "histhori_regular 2.0"
661!-
662  IF (rectilinear) THEN
663    ndim = 1
664    dims(1:1) = (/ xid(pfileid) /)
665  ENDIF
666  iret = NF90_DEF_VAR (ncid,lon_name,NF90_FLOAT,dims(1:ndim),nlonid)
667  iret = NF90_PUT_ATT (ncid,nlonid,'axis',"X")
668  iret = NF90_PUT_ATT (ncid,nlonid,'units',"degrees_east")
669  iret = NF90_PUT_ATT (ncid,nlonid,'valid_min', &
670 &                     REAL(MINVAL(plon),KIND=4))
671  iret = NF90_PUT_ATT (ncid,nlonid,'valid_max', &
672 &                     REAL(MAXVAL(plon),KIND=4))
673  iret = NF90_PUT_ATT (ncid,nlonid,'long_name',"Longitude")
674  iret = NF90_PUT_ATT (ncid,nlonid,'nav_model',TRIM(tmp_title))
675!-
676! 3.0 Latitude
677!-
678  IF (check) WRITE(*,*) "histhori_regular 3.0"
679!-
680  IF (rectilinear) THEN
681    ndim = 1
682    dims(1:1) = (/ yid(pfileid) /)
683  ENDIF
684  iret = NF90_DEF_VAR (ncid,lat_name,NF90_FLOAT,dims(1:ndim),nlatid)
685  iret = NF90_PUT_ATT (ncid,nlatid,'axis',"Y")
686  iret = NF90_PUT_ATT (ncid,nlatid,'units',"degrees_north")
687  iret = NF90_PUT_ATT (ncid,nlatid,'valid_min', &
688 &                     REAL(MINVAL(plat),KIND=4))
689  iret = NF90_PUT_ATT (ncid,nlatid,'valid_max', &
690 &                     REAL(MAXVAL(plat),KIND=4))
691  iret = NF90_PUT_ATT (ncid,nlatid,'long_name',"Latitude")
692  iret = NF90_PUT_ATT (ncid,nlatid,'nav_model',TRIM(tmp_title))
693!-
694  iret = NF90_ENDDEF (ncid)
695!-
696! 4.0 storing the geographical coordinates
697!-
698  IF (check) WRITE(*,*) "histhori_regular 4.0"
699!-
700  orix = slab_ori(pfileid,1)
701  oriy = slab_ori(pfileid,2)
702  par_szx = slab_sz(pfileid,1)
703  par_szy = slab_sz(pfileid,2)
704!-
705! Transfer the longitude
706!-
707  IF (rectilinear) THEN
708    iret = NF90_PUT_VAR (ncid,nlonid,plon(orix:orix+par_szx-1,1))
709  ELSE
710    iret = NF90_PUT_VAR (ncid,nlonid, &
711 &           plon(orix:orix+par_szx-1,oriy:oriy+par_szy-1))
712  ENDIF
713!-
714! Transfer the latitude
715!-
716  IF (rectilinear) THEN
717    iret = NF90_PUT_VAR (ncid,nlatid,plat(1,oriy:oriy+par_szy-1))
718  ELSE
719    iret = NF90_PUT_VAR (ncid,nlatid, &
720 &           plat(orix:orix+par_szx-1,oriy:oriy+par_szy-1))
721  ENDIF
722!-
723  iret = NF90_REDEF (ncid)
724!------------------------------
725END SUBROUTINE histhori_regular
726!===
727SUBROUTINE histhori_irregular &
728 &  (pfileid, pim, plon, plon_bounds, plat, plat_bounds, &
729 &   phname, phtitle, phid)
730!---------------------------------------------------------------------
731!- This subroutine is made to declare a new horizontale grid.
732!- It has to have the same number of points as
733!- the original and thus in this routine we will only
734!- add two variable (longitude and latitude).
735!- Any variable in the file can thus point to this pair
736!- through an attribute. This routine is very usefull
737!- to allow staggered grids.
738!-
739!- INPUT
740!-
741!- pfileid     : The id of the file to which the grid should be added
742!- pim         : Size in the longitude direction
743!- plon        : The longitudes
744!- plon_bounds : The boundaries of the grid in longitude
745!- plat        : The latitudes
746!- plat_bounds : Boundaries of the grid in latitude
747!- phname      : The name of grid
748!- phtitle     : The title of the grid
749!-
750!- OUTPUT
751!-
752!- phid    : Id of the created grid
753!---------------------------------------------------------------------
754  IMPLICIT NONE
755!-
756  INTEGER,INTENT(IN) :: pfileid, pim
757  REAL,DIMENSION(pim),INTENT(IN) :: plon, plat
758  REAL,DIMENSION(:,:),INTENT(IN) :: plon_bounds, plat_bounds
759  CHARACTER(LEN=*), INTENT(IN) :: phname, phtitle
760  INTEGER,INTENT(OUT) :: phid
761!-
762  CHARACTER(LEN=25) :: lon_name, lat_name
763  CHARACTER(LEN=30) :: lonbound_name, latbound_name
764  CHARACTER(LEN=80) :: tmp_title, tmp_name, longname
765  INTEGER :: ndim, dims(2)
766  INTEGER :: ndimb, dimsb(2)
767  INTEGER :: nbbounds
768  INTEGER :: nlonid, nlatid, nlonidb, nlatidb
769  INTEGER :: iret, ncid, twoid
770  LOGICAL :: transp = .FALSE.
771  REAL, ALLOCATABLE, DIMENSION(:,:) :: bounds_trans
772!-
773  LOGICAL :: check = .FALSE.
774!---------------------------------------------------------------------
775!-
776! 1.0 Check that all fits in the buffers
777!-
778  IF (    (pim /= full_size(pfileid,1)) &
779    & .OR.(full_size(pfileid,2) /= 1) ) THEN
780    CALL ipslerr (3,"histhori", &
781   &  'The new horizontal grid does not have the same size', &
782   &  'as the one provided to histbeg. This is not yet ', &
783   &  'possible in the hist package.')
784  ENDIF
785!-
786! 1.1 Create all the variables needed
787!-
788  IF (check) WRITE(*,*) 'histhori_irregular 1.0'
789!-
790  ncid = ncdf_ids(pfileid)
791!-
792  IF     (SIZE(plon_bounds,DIM=1) == pim) THEN
793    nbbounds = SIZE(plon_bounds,DIM=2)
794    transp = .TRUE.
795  ELSEIF (SIZE(plon_bounds,DIM=2) == pim) THEN
796    nbbounds = SIZE(plon_bounds,DIM=1)
797    transp = .FALSE.
798  ELSE
799    CALL ipslerr (3,"histhori", &
800   &  'The boundary variable does not have any axis corresponding', &
801   &  'to the size of the longitude or latitude variable', &
802   &  '.')
803  ENDIF
804!-
805  IF (.NOT.ALLOCATED(bounds_trans)) THEN
806    ALLOCATE(bounds_trans(nbbounds,pim))
807  ENDIF
808!-
809  iret = NF90_DEF_DIM (ncid, 'nbnd', nbbounds, twoid)
810  ndim = 1
811  dims(1) = xid(pfileid)
812  ndimb = 2
813  dimsb(1:2) = (/ twoid, xid(pfileid) /)
814!-
815  tmp_name = phname
816  IF (nb_hax(pfileid) == 0) THEN
817    lon_name = 'nav_lon'
818    lat_name = 'nav_lat'
819  ELSE
820    lon_name = 'nav_lon_'//TRIM(tmp_name)
821    lat_name = 'nav_lat_'//TRIM(tmp_name)
822  ENDIF
823  lonbound_name = TRIM(lon_name)//'_bounds'
824  latbound_name = TRIM(lat_name)//'_bounds'
825!-
826! 1.2 Save the informations
827!-
828  phid =  nb_hax(pfileid)+1
829  nb_hax(pfileid) = phid
830!-
831  hax_name(pfileid,phid,1:2) = (/ lon_name, lat_name /)
832  tmp_title = phtitle
833!-
834! 2.0 Longitude
835!-
836  IF (check) WRITE(*,*) "histhori_irregular 2.0"
837!-
838  iret = NF90_DEF_VAR (ncid,lon_name,NF90_FLOAT,dims(1:ndim),nlonid)
839  iret = NF90_PUT_ATT (ncid,nlonid,'axis',"X")
840  iret = NF90_PUT_ATT (ncid,nlonid,'units',"degrees_east")
841  iret = NF90_PUT_ATT (ncid,nlonid,'valid_min', &
842 &                     REAL(MINVAL(plon),KIND=4))
843  iret = NF90_PUT_ATT (ncid,nlonid,'valid_max', &
844 &                     REAL(MAXVAL(plon),KIND=4))
845  iret = NF90_PUT_ATT (ncid,nlonid,'long_name',"Longitude")
846  iret = NF90_PUT_ATT (ncid,nlonid,'nav_model',TRIM(tmp_title))
847!-
848! 2.1 Longitude bounds
849!-
850  iret = NF90_PUT_ATT (ncid,nlonid,'bounds',TRIM(lonbound_name))
851  iret = NF90_DEF_VAR (ncid,lonbound_name,NF90_FLOAT, &
852 &                     dimsb(1:ndimb),nlonidb)
853  longname = 'Boundaries for coordinate variable '//TRIM(lon_name)
854  iret = NF90_PUT_ATT (ncid,nlonidb,'long_name',TRIM(longname))
855!-
856! 3.0 Latitude
857!-
858  IF (check) WRITE(*,*) "histhori_irregular 3.0"
859!-
860  iret = NF90_DEF_VAR (ncid,lat_name,NF90_FLOAT,dims(1:ndim),nlatid)
861  iret = NF90_PUT_ATT (ncid,nlatid,'axis',"Y")
862  iret = NF90_PUT_ATT (ncid,nlatid,'units',"degrees_north")
863  iret = NF90_PUT_ATT (ncid,nlatid,'valid_min', &
864 &                     REAL(MINVAL(plat),KIND=4))
865  iret = NF90_PUT_ATT (ncid,nlatid,'valid_max', &
866 &                     REAL(MAXVAL(plat),KIND=4))
867  iret = NF90_PUT_ATT (ncid,nlatid,'long_name',"Latitude")
868  iret = NF90_PUT_ATT (ncid,nlatid,'nav_model',TRIM(tmp_title))
869!-
870! 3.1 Latitude bounds
871!-
872  iret = NF90_PUT_ATT (ncid,nlatid,'bounds',TRIM(latbound_name))
873  iret = NF90_DEF_VAR (ncid,latbound_name,NF90_FLOAT, &
874 &                     dimsb(1:ndimb),nlatidb)
875  longname = 'Boundaries for coordinate variable '//TRIM(lat_name)
876  iret = NF90_PUT_ATT (ncid,nlatidb,'long_name',TRIM(longname))
877!-
878  iret = NF90_ENDDEF (ncid)
879!-
880! 4.0 storing the geographical coordinates
881!-
882  IF (check) WRITE(*,*) "histhori_irregular 4.0"
883!-
884! 4.1 Write the longitude
885!-
886  iret = NF90_PUT_VAR (ncid, nlonid, plon(1:pim))
887!-
888! 4.2 Write the longitude bounds
889!-
890  IF (transp) THEN
891    bounds_trans = TRANSPOSE(plon_bounds)
892  ELSE
893    bounds_trans = plon_bounds
894  ENDIF
895  iret = NF90_PUT_VAR (ncid, nlonidb, bounds_trans(1:nbbounds,1:pim))
896!-
897! 4.3 Write the latitude
898!-
899  iret = NF90_PUT_VAR (ncid, nlatid, plat(1:pim))
900!-
901! 4.4 Write the latitude bounds
902!-
903  IF (transp) THEN
904    bounds_trans = TRANSPOSE(plat_bounds)
905  ELSE
906    bounds_trans = plat_bounds
907  ENDIF
908  iret = NF90_PUT_VAR (ncid,nlatidb,bounds_trans(1:nbbounds,1:pim))
909!-
910  iret = NF90_REDEF (ncid)
911!--------------------------------
912END SUBROUTINE histhori_irregular
913!===
914SUBROUTINE histvert (pfileid,pzaxname,pzaxtitle,pzaxunit, &
915 &                   pzsize,pzvalues,pzaxid,pdirect)
916!---------------------------------------------------------------------
917!- This subroutine defines a vertical axis and returns it s id.
918!- It gives the user the possibility to the user to define many
919!- different vertical axes. For each variable defined with histdef a
920!- vertical axis can be specified with by it s ID.
921!-
922!- INPUT
923!-
924!- pfileid  : ID of the file the variable should be archived in
925!- pzaxname : Name of the vertical axis
926!- pzaxtitle: title of the vertical axis
927!- pzaxunit : Units of the vertical axis (no units ih blank string)
928!- pzsize   : size of the vertical axis
929!- pzvalues : Coordinate values of the vetical axis
930!-
931!- pdirect  : is an optional argument which allows to specify the
932!-            the positive direction of the axis : up or down.
933!- OUTPUT
934!-
935!- pzaxid   : Returns the ID of the axis.
936!-            Note that this is not the netCDF ID !
937!-
938!- VERSION
939!-
940!---------------------------------------------------------------------
941  IMPLICIT NONE
942!-
943  INTEGER,INTENT(IN) :: pfileid,pzsize
944  CHARACTER(LEN=*),INTENT(IN) :: pzaxname,pzaxunit,pzaxtitle
945  REAL,INTENT(IN) :: pzvalues(pzsize)
946  INTEGER, INTENT(OUT) :: pzaxid
947  CHARACTER(LEN=*),INTENT(IN),OPTIONAL :: pdirect
948!-
949  INTEGER :: pos,iv,zdimid,zaxid_tmp
950  CHARACTER(LEN=70) :: str71
951  CHARACTER(LEN=80) :: str80
952  CHARACTER(LEN=20) :: direction
953  INTEGER :: iret, leng, ncid
954  LOGICAL :: check = .FALSE.
955!---------------------------------------------------------------------
956!-
957! 1.0 Verifications :
958!    Do we have enough space for an extra axis ?
959!    Is the name already in use ?
960!-
961  IF (check) WRITE(*,*) "histvert : 1.0 Verifications", &
962 &                      pzaxname,'---',pzaxunit,'---',pzaxtitle
963!-
964! - Direction of axis. Can we get if from the user.
965!   If not we put unknown.
966!-
967  IF (PRESENT(pdirect)) THEN
968    direction = TRIM(pdirect)
969    CALL strlowercase (direction)
970  ELSE
971    direction = 'unknown'
972  ENDIF
973!-
974! Check the consistency of the attribute
975!-
976  IF (     (direction /= 'unknown') &
977 &    .AND.(direction /= 'up')      &
978 &    .AND.(direction /= 'down')   ) THEN
979    direction = 'unknown'
980    str80 = 'The specified axis was : '//TRIM(direction)
981    CALL ipslerr (2,"histvert",&
982   & "The specified direction for the vertical axis is not possible.",&
983   & "it is replaced by : unknown", str80)
984  ENDIF
985!-
986  IF (nb_zax(pfileid)+1 > nb_zax_max) THEN
987    CALL ipslerr (3,"histvert", &
988   &  'Table of vertical axes too small. You should increase ',&
989   &  'nb_zax_max in histcom.f90 in order to accomodate all ', &
990   &  'these variables ')
991  ENDIF
992!-
993  iv = nb_zax(pfileid)
994  IF (iv > 1) THEN
995    CALL find_str (zax_name(pfileid,1:iv-1),pzaxname,pos)
996  ELSE
997    pos = 0
998  ENDIF
999!-
1000  IF (pos > 0) THEN
1001    WRITE(str71,'("Check variable ",A," in file",I3)') &
1002 &    TRIM(pzaxname),pfileid
1003    CALL ipslerr (3,"histvert", &
1004 &    "Vertical axis already exists",TRIM(str71), &
1005 &    "Can also be a wrong file ID in another declaration")
1006  ENDIF
1007!-
1008  iv = nb_zax(pfileid)+1
1009!-
1010! 2.0 Add the information to the file
1011!-
1012  IF (check) &
1013 &  WRITE(*,*) "histvert : 2.0 Add the information to the file"
1014!-
1015  ncid = ncdf_ids(pfileid)
1016!-
1017  leng = MIN(LEN_TRIM(pzaxname),20)
1018  iret = NF90_DEF_DIM (ncid,pzaxname(1:leng),pzsize,zaxid_tmp)
1019  iret = NF90_DEF_VAR (ncid,pzaxname(1:leng),NF90_FLOAT, &
1020 &                     zaxid_tmp,zdimid)
1021  iret = NF90_PUT_ATT (ncid,zdimid,'axis',"Z")
1022  leng = MIN(LEN_TRIM(pzaxunit),20)
1023  IF (leng > 0) THEN
1024    iret = NF90_PUT_ATT (ncid,zdimid,'units',pzaxunit(1:leng))
1025  ENDIF
1026  iret = NF90_PUT_ATT (ncid,zdimid,'positive',TRIM(direction))
1027  iret = NF90_PUT_ATT (ncid,zdimid,'valid_min', &
1028 &                     REAL(MINVAL(pzvalues(1:pzsize)),KIND=4))
1029  iret = NF90_PUT_ATT (ncid,zdimid,'valid_max', &
1030 &                     REAL(MAXVAL(pzvalues(1:pzsize)),KIND=4))
1031  leng = MIN(LEN_TRIM(pzaxname),20)
1032  iret = NF90_PUT_ATT (ncid,zdimid,'title',pzaxname(1:leng))
1033  leng = MIN(LEN_TRIM(pzaxtitle),80)
1034  iret = NF90_PUT_ATT (ncid,zdimid,'long_name',pzaxtitle(1:leng))
1035!-
1036  iret = NF90_ENDDEF (ncid)
1037!-
1038  iret = NF90_PUT_VAR (ncid,zdimid,pzvalues(1:pzsize))
1039!-
1040  iret = NF90_REDEF (ncid)
1041!-
1042!- 3.0 add the information to the common
1043!-
1044  IF (check) &
1045  &  WRITE(*,*) "histvert : 3.0 add the information to the common"
1046!-
1047  nb_zax(pfileid) = iv
1048  zax_size(pfileid, iv) = pzsize
1049  zax_name(pfileid, iv) = pzaxname
1050  zax_ids(pfileid, iv) = zaxid_tmp
1051  pzaxid =  iv
1052!----------------------
1053END SUBROUTINE histvert
1054!===
1055SUBROUTINE histdef (pfileid, pvarname, ptitle, punit, &
1056 &                  pxsize, pysize, phoriid, pzsize,  &
1057 &                  par_oriz, par_szz, pzid,          &
1058 &                  pnbbyt, popp, pfreq_opp, pfreq_wrt, var_range)
1059!---------------------------------------------------------------------
1060!- With this subroutine each variable to be archived on the history
1061!- tape should be declared.
1062!-
1063!- It gives the user the choise of operation
1064!- to be performed on the variables, the frequency of this operation
1065!- and finaly the frequency of the archiving.
1066!-
1067!- INPUT
1068!-
1069!- pfileid  : ID of the file the variable should be archived in
1070!- pvarname : Name of the variable, short and easy to remember
1071!- ptitle   : Full name of the variable
1072!- punit    : Units of the variable (no units if blank string)
1073!-
1074!- The next 3 arguments give the size of that data
1075!- that will be passed to histwrite. The zoom will be
1076!- done there with the horizontal information obtained
1077!- in histbeg and the vertical information to follow.
1078!-
1079!- pxsize   : Size in X direction (size of the data that will be
1080!-            given to histwrite)
1081!- pysize   : Size in Y direction
1082!- phoriid  : ID of the horizontal axis
1083!-
1084!- The next two arguments give the vertical zoom to use.
1085!-
1086!- pzsize   : Size in Z direction (If 1 then no axis is declared
1087!-            for this variable and pzid is not used)
1088!- par_oriz : Off set of the zoom
1089!- par_szz  : Size of the zoom
1090!-
1091!- pzid     : ID of the vertical axis to use. It has to have
1092!-            the size of the zoom.
1093!- pnbbyt   : Number of bytes on which to store in netCDF (Not opp.)
1094!- popp     : Operation to be performed. The following options
1095!-            exist today :
1096!- inst : keeps instantaneous values for writting
1097!- ave  : Computes the average from call between writes
1098!- pfreq_opp: Frequency of this operation (in seconds)
1099!- pfreq_wrt: Frequency at which the variable should be
1100!-            written (in seconds)
1101!- var_range: Range of the variable.
1102!-            If the minimum is greater than the maximum,
1103!-            the values will be calculated.
1104!-
1105!- VERSION
1106!---------------------------------------------------------------------
1107  IMPLICIT NONE
1108!-
1109  INTEGER,INTENT(IN) :: pfileid, pxsize, pysize, pzsize, pzid
1110  INTEGER,INTENT(IN) :: par_oriz, par_szz, pnbbyt, phoriid
1111  CHARACTER(LEN=*),INTENT(IN) :: pvarname, punit, popp
1112  CHARACTER(LEN=*),INTENT(IN) :: ptitle
1113  REAL,INTENT(IN) :: pfreq_opp, pfreq_wrt
1114  REAL,DIMENSION(2),OPTIONAL,INTENT(IN) :: var_range
1115!-
1116  INTEGER :: iv,i
1117  CHARACTER(LEN=70) :: str70, str71, str72
1118  CHARACTER(LEN=20) :: tmp_name
1119  CHARACTER(LEN=40) :: str40
1120  CHARACTER(LEN=10) :: str10
1121  CHARACTER(LEN=80) :: tmp_str80
1122  CHARACTER(LEN=7) :: tmp_topp, tmp_sopp(nbopp_max)
1123  CHARACTER(LEN=120) :: ex_topps
1124  REAL :: tmp_scal(nbopp_max), un_an, un_jour, test_fopp, test_fwrt
1125  INTEGER :: pos, buff_sz
1126!-
1127  LOGICAL :: check = .FALSE.
1128!---------------------------------------------------------------------
1129  ex_topps = 'ave, inst, t_min, t_max, t_sum, once, never, l_max, l_min'
1130!-
1131  nb_var(pfileid) = nb_var(pfileid)+1
1132  iv = nb_var(pfileid)
1133!-
1134  IF (iv > nb_var_max) THEN
1135    CALL ipslerr (3,"histdef", &
1136   &  'Table of variables too small. You should increase nb_var_max',&
1137   &  'in histcom.f90 in order to accomodate all these variables', &
1138   &  ' ')
1139  ENDIF
1140!-
1141! 1.0 Transfer informations on the variable to the common
1142!     and verify that it does not already exist
1143!-
1144  IF (check) WRITE(*,*) "histdef : 1.0"
1145!-
1146  IF (iv > 1) THEN
1147    CALL find_str (name(pfileid,1:iv-1),pvarname,pos)
1148  ELSE
1149    pos = 0
1150  ENDIF
1151!-
1152  IF (pos > 0) THEN
1153    str70 = "Variable already exists"
1154    WRITE(str71,'("Check variable  ",a," in file",I3)') &
1155 &    TRIM(pvarname),pfileid
1156    str72 = "Can also be a wrong file ID in another declaration"
1157    CALL ipslerr (3,"histdef", str70, str71, str72)
1158  ENDIF
1159!-
1160  name(pfileid,iv) = pvarname
1161  title(pfileid,iv) = ptitle
1162  unit_name(pfileid,iv) = punit
1163  tmp_name =  name(pfileid,iv)
1164!-
1165! 1.1 decode the operations
1166!-
1167  fullop(pfileid,iv) = popp
1168  tmp_str80 = popp
1169  CALL buildop &
1170 &  (tmp_str80, ex_topps, tmp_topp, nbopp_max, missing_val, &
1171 &   tmp_sopp, tmp_scal, nbopp(pfileid,iv))
1172!-
1173  topp(pfileid,iv) = tmp_topp
1174  DO i=1,nbopp(pfileid,iv)
1175    sopps(pfileid,iv,i) = tmp_sopp(i)
1176    scal(pfileid,iv,i) = tmp_scal(i)
1177  ENDDO
1178!-
1179! 1.2 If we have an even number of operations
1180!     then we need to add identity
1181!-
1182  IF (2*INT(nbopp(pfileid,iv)/2.0) == nbopp(pfileid,iv)) THEN
1183    nbopp(pfileid,iv) = nbopp(pfileid,iv)+1
1184    sopps(pfileid,iv,nbopp(pfileid,iv)) = 'ident'
1185    scal(pfileid,iv,nbopp(pfileid,iv)) = missing_val
1186  ENDIF
1187!-
1188! 2.0 Put the size of the variable in the common and check
1189!-
1190  IF (check) &
1191 &  WRITE(*,*) "histdef : 2.0", pfileid,iv,nbopp(pfileid,iv), &
1192 &    sopps(pfileid,iv,1:nbopp(pfileid,iv)), &
1193 &    scal(pfileid,iv,1:nbopp(pfileid,iv))
1194!-
1195  scsize(pfileid,iv,1:3) = (/ pxsize, pysize, pzsize /)
1196!-
1197  zorig(pfileid,iv,1:3) = &
1198 &  (/ slab_ori(pfileid,1), slab_ori(pfileid,2), par_oriz /)
1199!-
1200  zsize(pfileid,iv,1:3) = &
1201 &  (/ slab_sz(pfileid,1), slab_sz(pfileid,2), par_szz /)
1202!-
1203! Is the size of the full array the same as that of the coordinates  ?
1204!-
1205  IF (    (pxsize > full_size(pfileid,1)) &
1206 &    .OR.(pysize > full_size(pfileid,2)) ) THEN
1207!-
1208    str70 = "The size of the variable is different "// &
1209 &          "from the one of the coordinates"
1210    WRITE(str71,'("Size of coordinates :", 2I4)') &
1211 &   full_size(pfileid,1), full_size(pfileid,2)
1212    WRITE(str72,'("Size declared for variable ",a," :",2I4)') &
1213 &   TRIM(tmp_name), pxsize, pysize
1214    CALL ipslerr (3,"histdef", str70, str71, str72)
1215  ENDIF
1216!-
1217! Is the size of the zoom smaler than the coordinates ?
1218!-
1219  IF (    (full_size(pfileid,1) < slab_sz(pfileid,1)) &
1220 &    .OR.(full_size(pfileid,2) < slab_sz(pfileid,2)) ) THEN
1221    str70 = &
1222 &   "Size of variable should be greater or equal to those of the zoom"
1223    WRITE(str71,'("Size of XY zoom :", 2I4)') &
1224 &   slab_sz(pfileid,1),slab_sz(pfileid,1)
1225    WRITE(str72,'("Size declared for variable ",a," :",2I4)') &
1226 &   TRIM(tmp_name), pxsize, pysize
1227    CALL ipslerr (3,"histdef", str70, str71, str72)
1228  ENDIF
1229!-
1230! 2.1 We store the horizontal grid information with minimal
1231!     and a fall back onto the default grid
1232!-
1233  IF ( (phoriid > 0).AND.(phoriid <= nb_hax(pfileid)) ) THEN
1234    var_haxid(pfileid,iv) = phoriid
1235  ELSE
1236    var_haxid(pfileid,iv) = 1
1237    CALL ipslerr (2,"histdef", &
1238   &  'We use the default grid for variable as an invalide',&
1239   &  'ID was provided for variable : ', pvarname)
1240  ENDIF
1241!-
1242! 2.2 Check the vertical coordinates if needed
1243!-
1244  IF (par_szz > 1) THEN
1245!-
1246!-- Does the vertical coordinate exist ?
1247!-
1248    IF (pzid > nb_zax(pfileid)) THEN
1249      WRITE(str70, &
1250 &    '("The vertical coordinate chosen for variable ",a)') &
1251 &     TRIM(tmp_name)
1252      str71 = " Does not exist."
1253      CALL ipslerr (3,"histdef",str70,str71, " ")
1254    ENDIF
1255!-
1256!-- Is the vertical size of the variable equal to that of the axis ?
1257!-
1258    IF (par_szz /= zax_size(pfileid,pzid)) THEN
1259      str70 = "The size of the zoom does not correspond "// &
1260 &            "to the size of the chosen vertical axis"
1261      WRITE(str71,'("Size of zoom in z :", I4)') par_szz
1262      WRITE(str72,'("Size declared for axis ",A," :",I4)') &
1263 &     TRIM(zax_name(pfileid,pzid)),zax_size(pfileid,pzid)
1264      CALL ipslerr (3,"histdef", str70, str71, str72)
1265    ENDIF
1266!-
1267!-- Is the zoom smaler that the total size of the variable ?
1268!-
1269    IF (pzsize < par_szz) THEN
1270      str70 = "The vertical size of variable "// &
1271 &            "is smaller than that of the zoom."
1272      WRITE(str71,'("Declared vertical size of data :", I5)') pzsize
1273      WRITE(str72,'("Size of zoom for variable ",a," = ",I5)') &
1274 &     TRIM(tmp_name),par_szz
1275      CALL ipslerr (3,"histdef", str70, str71, str72)
1276    ENDIF
1277    var_zaxid(pfileid,iv) = pzid
1278  ELSE
1279    var_zaxid(pfileid,iv) = -99
1280  ENDIF
1281!-
1282! 3.0 Determine the position of the variable in the buffer
1283!     If it is instantaneous output then we do not use the buffer
1284!-
1285  IF (check) WRITE(*,*) "histdef : 3.0"
1286!-
1287! 3.1 We get the size of the arrays histwrite will get and check
1288!     that they fit into the tmp_buffer
1289!-
1290  buff_sz = zsize(pfileid,iv,1)*zsize(pfileid,iv,2)*zsize(pfileid,iv,3)
1291!-
1292! 3.2 move the pointer of the buffer array for operation
1293!     which need bufferisation
1294!-
1295  IF (     (TRIM(tmp_topp) /= "inst") &
1296 &    .AND.(TRIM(tmp_topp) /= "once") &
1297 &    .AND.(TRIM(tmp_topp) /= "never") )THEN
1298    point(pfileid,iv) = buff_pos+1
1299    buff_pos = buff_pos+buff_sz
1300    IF (check) THEN
1301      WRITE(*,*) "histdef : 3.2 bufpos for iv = ",iv, &
1302 &               " pfileid = ",pfileid," is = ",point(pfileid,iv)
1303    ENDIF
1304  ENDIF
1305!-
1306! 4.0 Transfer the frequency of the operations and check
1307!     for validity. We have to pay attention to negative values
1308!     of the frequency which indicate monthly time-steps.
1309!     The strategy is to bring it back to seconds for the tests
1310!-
1311  IF (check) WRITE(*,*) "histdef : 4.0"
1312!-
1313  freq_opp(pfileid,iv) = pfreq_opp
1314  freq_wrt(pfileid,iv) = pfreq_wrt
1315!-
1316  CALL ioget_calendar(un_an, un_jour)
1317  IF (pfreq_opp < 0) THEN
1318    CALL ioget_calendar(un_an)
1319    test_fopp = pfreq_opp*(-1.)*un_an/12.*un_jour
1320  ELSE
1321    test_fopp = pfreq_opp
1322  ENDIF
1323  IF (pfreq_wrt < 0) THEN
1324    CALL ioget_calendar(un_an)
1325    test_fwrt = pfreq_wrt*(-1.)*un_an/12.*un_jour
1326  ELSE
1327    test_fwrt = pfreq_wrt
1328  ENDIF
1329!-
1330! 4.1 Frequency of operations and output should be larger than deltat !
1331!-
1332  IF (test_fopp < deltat(pfileid)) THEN
1333    str70 = 'Frequency of operations should be larger than deltat'
1334    WRITE(str71,'("It is not the case for variable ",a," :",F10.4)') &
1335 &   TRIM(tmp_name),pfreq_opp
1336    str72 = "PATCH : frequency set to deltat"
1337!-
1338    CALL ipslerr (2,"histdef", str70, str71, str72)
1339!-
1340    freq_opp(pfileid,iv) = deltat(pfileid)
1341  ENDIF
1342!-
1343  IF (test_fwrt < deltat(pfileid)) THEN
1344    str70 = 'Frequency of output should be larger than deltat'
1345    WRITE(str71,'("It is not the case for variable ",a," :",F10.4)') &
1346 &   TRIM(tmp_name),pfreq_wrt
1347    str72 = "PATCH : frequency set to deltat"
1348!-
1349    CALL ipslerr (2,"histdef", str70, str71, str72)
1350!-
1351    freq_wrt(pfileid,iv) = deltat(pfileid)
1352  ENDIF
1353!-
1354! 4.2 First the existence of the operation is tested and then
1355!     its compaticility with the choice of frequencies
1356!-
1357  IF (TRIM(tmp_topp) == "inst") THEN
1358    IF (test_fopp /= test_fwrt) THEN
1359      str70 = 'For instantaneous output the frequency '// &
1360 &            'of operations and output'
1361      WRITE(str71, &
1362 &     '("should be the same, this was not case for variable ",a)') &
1363 &      TRIM(tmp_name)
1364      str72 = "PATCH : The smalest frequency of both is used"
1365      CALL ipslerr (2,"histdef", str70, str71, str72)
1366      IF (test_fopp < test_fwrt) THEN
1367        freq_opp(pfileid,iv) = pfreq_opp
1368        freq_wrt(pfileid,iv) = pfreq_opp
1369      ELSE
1370        freq_opp(pfileid,iv) = pfreq_wrt
1371        freq_wrt(pfileid,iv) = pfreq_wrt
1372      ENDIF
1373    ENDIF
1374  ELSE IF (INDEX(ex_topps,TRIM(tmp_topp)) > 0) THEN
1375    IF (test_fopp > test_fwrt) THEN
1376      str70 = 'For averages the frequency of operations '// &
1377&             'should be smaller or equal'
1378      WRITE(str71, &
1379 &     '("to that of output. It is not the case for variable ",a)') &
1380 &     TRIM(tmp_name)
1381      str72 = 'PATCH : The output frequency is used for both'
1382      CALL ipslerr (2,"histdef", str70, str71, str72)
1383      freq_opp(pfileid,iv) = pfreq_wrt
1384    ENDIF
1385  ELSE
1386    WRITE (str70,'("Operation on variable ",a," is unknown")') &
1387&    TRIM(tmp_name)
1388    WRITE (str71, '("operation requested is :",a)') tmp_topp
1389    WRITE (str72, '("File ID :",I3)') pfileid
1390    CALL ipslerr (3,"histdef", str70, str71, str72)
1391  ENDIF
1392!-
1393! 5.0 Initialize other variables of the common
1394!-
1395  IF (check) WRITE(*,*) "histdef : 5.0"
1396!-
1397  hist_wrt_rng(pfileid,iv) = (PRESENT(var_range))
1398  IF (hist_wrt_rng(pfileid,iv)) THEN
1399    hist_calc_rng(pfileid,iv) = (var_range(1) > var_range(2))
1400    IF (hist_calc_rng(pfileid,iv)) THEN
1401      hist_minmax(pfileid,iv,1:2) = &
1402 &      (/ ABS(missing_val), -ABS(missing_val) /)
1403    ELSE
1404      hist_minmax(pfileid,iv,1:2) = var_range(1:2)
1405    ENDIF
1406  ENDIF
1407!-
1408! - freq_opp(pfileid,iv)/2./deltat(pfileid)
1409  last_opp(pfileid,iv) = itau0(pfileid)
1410! - freq_wrt(pfileid,iv)/2./deltat(pfileid)
1411  last_wrt(pfileid,iv) = itau0(pfileid)
1412! - freq_opp(pfileid,iv)/2./deltat(pfileid)
1413  last_opp_chk(pfileid,iv) = itau0(pfileid)
1414! - freq_wrt(pfileid,iv)/2./deltat(pfileid)
1415  last_wrt_chk(pfileid,iv) = itau0(pfileid)
1416  nb_opp(pfileid,iv) = 0
1417  nb_wrt(pfileid,iv) = 0
1418!-
1419! 6.0 Get the time axis for this variable
1420!-
1421  IF (check) WRITE(*,*) "histdef : 6.0"
1422!-
1423  IF (freq_wrt(pfileid,iv) > 0) THEN
1424    WRITE(str10,'(I8.8)') INT(freq_wrt(pfileid,iv))
1425    str40 = TRIM(tmp_topp)//"_"//TRIM(str10)
1426  ELSE
1427    WRITE(str10,'(I2.2,"month")') ABS(INT(freq_wrt(pfileid,iv)))
1428    str40 = TRIM(tmp_topp)//"_"//TRIM(str10)
1429  ENDIF
1430  CALL find_str (tax_name(pfileid,1:nb_tax(pfileid)),str40,pos)
1431!-
1432! No time axis for once, l_max, l_min or never operation
1433!-
1434  IF (     (TRIM(tmp_topp) /= 'once')  &
1435 &    .AND.(TRIM(tmp_topp) /= 'never') &
1436 &    .AND.(TRIM(tmp_topp) /= 'l_max') &
1437 &    .AND.(TRIM(tmp_topp) /= 'l_min') ) THEN
1438    IF (pos < 0) THEN
1439      nb_tax(pfileid) = nb_tax(pfileid)+1
1440      tax_name(pfileid,nb_tax(pfileid)) = str40
1441      tax_last(pfileid,nb_tax(pfileid)) = 0
1442      var_axid(pfileid,iv) = nb_tax(pfileid)
1443    ELSE
1444      var_axid(pfileid,iv) = pos
1445    ENDIF
1446  ELSE
1447    IF (check)   WRITE(*,*) "histdef : 7.0 ",TRIM(tmp_topp),'----'
1448    var_axid(pfileid,iv) = -99
1449  ENDIF
1450!-
1451! 7.0 prepare frequence of writing and operation
1452!     for never or once operation
1453!-
1454  IF (    (TRIM(tmp_topp) == 'once')  &
1455 &    .OR.(TRIM(tmp_topp) == 'never') ) THEN
1456    freq_opp(pfileid,iv) = 0.
1457    freq_wrt(pfileid,iv) = 0.
1458  ENDIF
1459!---------------------
1460END SUBROUTINE histdef
1461!===
1462SUBROUTINE histend (pfileid)
1463!---------------------------------------------------------------------
1464!- This subroutine end the decalaration of variables and sets the
1465!- time axes in the netcdf file and puts it into the write mode.
1466!-
1467!- INPUT
1468!-
1469!- pfileid : ID of the file to be worked on
1470!-
1471!- VERSION
1472!-
1473!---------------------------------------------------------------------
1474  IMPLICIT NONE
1475!-
1476  INTEGER, INTENT(IN) :: pfileid
1477!-
1478  INTEGER :: ncid, ncvarid
1479  INTEGER :: iret, ndim, iv, itx, ziv
1480  INTEGER :: itax
1481  INTEGER :: dims(4), dim_cnt
1482  INTEGER :: year, month, day, hours, minutes
1483  REAL :: sec
1484  REAL :: rtime0
1485  CHARACTER(LEN=20) :: tname, tunit
1486  CHARACTER(LEN=30) :: str30
1487  CHARACTER(LEN=80) :: ttitle
1488  CHARACTER(LEN=120) :: assoc
1489  CHARACTER(LEN=70) :: str70
1490  CHARACTER(LEN=3),DIMENSION(12) :: cal =   &
1491 &  (/ 'JAN','FEB','MAR','APR','MAY','JUN', &
1492 &     'JUL','AUG','SEP','OCT','NOV','DEC' /)
1493  CHARACTER(LEN=7) :: tmp_opp
1494!-
1495  LOGICAL :: check = .FALSE.
1496!---------------------------------------------------------------------
1497  ncid = ncdf_ids(pfileid)
1498!-
1499! 1.0 Create the time axes
1500!-
1501  IF (check) WRITE(*,*) "histend : 1.0"
1502!---
1503  iret = NF90_DEF_DIM (ncid,'time_counter',NF90_UNLIMITED,tid(pfileid))
1504!-
1505! 1.1 Define all the time axes needed for this file
1506!-
1507  DO itx=1,nb_tax(pfileid)
1508    dims(1) = tid(pfileid)
1509    IF (nb_tax(pfileid) > 1) THEN
1510      str30 = "t_"//tax_name(pfileid,itx)
1511    ELSE
1512      str30 = "time_counter"
1513    ENDIF
1514    iret = NF90_DEF_VAR (ncid,str30,NF90_FLOAT, &
1515 &                       dims(1),tdimid(pfileid,itx))
1516    iret = NF90_PUT_ATT (ncid,tdimid(pfileid,itx),'axis',"T")
1517!---
1518!   To transform the current itau into a real date and take it
1519!   as the origin of the file requires the time counter to change.
1520!   Thus it is an operation the user has to ask for.
1521!   This function should thus only be re-instated
1522!   if there is a ioconf routine to control it.
1523!---
1524!-- rtime0 = itau2date(itau0(pfileid), date0(pfileid), deltat(pfileid))
1525    rtime0 = date0(pfileid)
1526!-
1527    CALL ju2ymds(rtime0, year, month, day, sec)
1528!---
1529!   Catch any error induced by a change in calendar !
1530!---
1531    IF (year < 0) THEN
1532      year = 2000+year
1533    ENDIF
1534!-
1535    hours = INT(sec/(60.*60.))
1536    minutes = INT((sec-hours*60.*60.)/60.)
1537    sec = sec-(hours*60.*60.+minutes*60.)
1538!-
1539    WRITE (UNIT=str70, &
1540 &   FMT='(A,I4.4,"-",I2.2,"-",I2.2," ",I2.2,":",I2.2,":",I2.2)') &
1541 &    'seconds since ',year,month,day,hours,minutes,INT(sec)
1542    iret = NF90_PUT_ATT (ncid,tdimid(pfileid,itx),'units',TRIM(str70))
1543!-
1544    CALL ioget_calendar (str30)
1545    iret = NF90_PUT_ATT (ncid,tdimid(pfileid,itx), &
1546 &                       'calendar',TRIM(str30))
1547!-
1548    iret = NF90_PUT_ATT (ncid,tdimid(pfileid,itx),'title','Time')
1549!-
1550    iret = NF90_PUT_ATT (ncid,tdimid(pfileid,itx), &
1551 &                       'long_name','Time axis')
1552!-
1553    WRITE (UNIT=str70, &
1554 &   FMT='(" ",I4.4,"-",A3,"-",I2.2," ",I2.2,":",I2.2,":",I2.2)') &
1555 &    year,cal(month),day,hours,minutes,INT(sec)
1556    iret = NF90_PUT_ATT (ncid,tdimid(pfileid,itx), &
1557 &                       'time_origin',TRIM(str70))
1558  ENDDO
1559!-
1560! 2.0 declare the variables
1561!-
1562  IF (check) WRITE(*,*) "histend : 2.0"
1563!-
1564  DO iv=1,nb_var(pfileid)
1565!---
1566    itax = var_axid(pfileid,iv)
1567!---
1568    tname = name(pfileid,iv)
1569    tunit = unit_name(pfileid,iv)
1570    ttitle = title(pfileid,iv)
1571!---
1572    IF (regular(pfileid) ) THEN
1573      dims(1:2) = (/ xid(pfileid), yid(pfileid) /)
1574      dim_cnt = 2
1575    ELSE
1576      dims(1) = xid(pfileid)
1577      dim_cnt = 1
1578    ENDIF
1579!---
1580    tmp_opp = topp(pfileid,iv)
1581    ziv = var_zaxid(pfileid,iv)
1582!---
1583!   2.1 dimension of field
1584!---
1585    IF ((TRIM(tmp_opp) /= 'never')) THEN
1586      IF (     (TRIM(tmp_opp) /= 'once')  &
1587     &    .AND.(TRIM(tmp_opp) /= 'l_max') &
1588     &    .AND.(TRIM(tmp_opp) /= 'l_min') ) THEN
1589        IF (ziv == -99) THEN
1590          ndim = dim_cnt+1
1591          dims(dim_cnt+1:dim_cnt+2) = (/ tid(pfileid), 0 /)
1592        ELSE
1593          ndim = dim_cnt+2
1594          dims(dim_cnt+1:dim_cnt+2) = (/ zax_ids(pfileid,ziv), tid(pfileid) /)
1595        ENDIF
1596      ELSE
1597        IF (ziv == -99) THEN
1598          ndim = dim_cnt
1599          dims(dim_cnt+1:dim_cnt+2) = (/ 0, 0 /)
1600        ELSE
1601          ndim = dim_cnt+1
1602          dims(dim_cnt+1:dim_cnt+2) = (/ zax_ids(pfileid,ziv), 0 /)
1603        ENDIF
1604      ENDIF
1605!-
1606      iret = NF90_DEF_VAR (ncid,TRIM(tname),NF90_FLOAT, &
1607 &                         dims(1:ABS(ndim)),ncvarid)
1608!-
1609      ncvar_ids(pfileid,iv) = ncvarid
1610!-
1611      IF (LEN_TRIM(tunit) > 0) THEN
1612        iret = NF90_PUT_ATT (ncid,ncvarid,'units',TRIM(tunit))
1613      ENDIF
1614!-
1615      iret = NF90_PUT_ATT (ncid,ncvarid,'missing_value', &
1616 &                         REAL(missing_val,KIND=4))
1617      IF (hist_wrt_rng(pfileid,iv)) THEN
1618        iret = NF90_PUT_ATT (ncid,ncvarid,'valid_min', &
1619 &                           REAL(hist_minmax(pfileid,iv,1),KIND=4))
1620        iret = NF90_PUT_ATT (ncid,ncvarid,'valid_max', &
1621 &                           REAL(hist_minmax(pfileid,iv,2),KIND=4))
1622      ENDIF
1623!-
1624      iret = NF90_PUT_ATT (ncid,ncvarid,'long_name',TRIM(ttitle))
1625!-
1626      iret = NF90_PUT_ATT (ncid,ncvarid,'short_name',TRIM(tname))
1627!-
1628      iret = NF90_PUT_ATT (ncid,ncvarid,'online_operation', &
1629 &                         TRIM(fullop(pfileid,iv)))
1630!-
1631      SELECT CASE(ndim)
1632      CASE(-3,2:4)
1633      CASE DEFAULT
1634        CALL ipslerr (3,"histend", &
1635       &  'less than 2 or more than 4 dimensions are not', &
1636       &  'allowed at this stage',' ')
1637      END SELECT
1638!-
1639      assoc='nav_lat nav_lon'
1640      ziv = var_zaxid(pfileid, iv)
1641      IF (ziv > 0) THEN
1642        str30 = zax_name(pfileid,ziv)
1643        assoc = TRIM(str30)//' '//TRIM(assoc)
1644      ENDIF
1645!-
1646      IF (itax > 0) THEN
1647        IF (nb_tax(pfileid) > 1) THEN
1648          str30 = "t_"//tax_name(pfileid,itax)
1649        ELSE
1650          str30 = "time_counter"
1651        ENDIF
1652        assoc = TRIM(str30)//' '//TRIM(assoc)
1653!-
1654        IF (check) THEN
1655          WRITE(*,*) "histend : 2.0.n, freq_opp, freq_wrt", &
1656 &                   freq_opp(pfileid,iv), freq_wrt(pfileid,iv)
1657        ENDIF
1658!-
1659        iret = NF90_PUT_ATT (ncid,ncvarid,'interval_operation', &
1660 &                           REAL(freq_opp(pfileid,iv),KIND=4))
1661        iret = NF90_PUT_ATT (ncid,ncvarid,'interval_write', &
1662 &                           REAL(freq_wrt(pfileid,iv),KIND=4))
1663      ENDIF
1664      iret = NF90_PUT_ATT (ncid,ncvarid,'coordinates',TRIM(assoc))
1665    ENDIF
1666  ENDDO
1667!-
1668! 2.2 Add DOMAIN attributes if needed
1669!-
1670  IF (dom_id_svg(pfileid) >= 0) THEN
1671    CALL flio_dom_att (ncid,dom_id_svg(pfileid))
1672  ENDIF
1673!-
1674! 3.0 Put the netcdf file into write mode
1675!-
1676  IF (check) WRITE(*,*) "histend : 3.0"
1677!-
1678  iret = NF90_ENDDEF (ncid)
1679!-
1680! 4.0 Give some informations to the user
1681!-
1682  IF (check) WRITE(*,*) "histend : 4.0"
1683!-
1684  WRITE(str70,'("All variables have been initialized on file :",I3)') pfileid
1685  CALL ipslerr (1,'histend',str70,'',' ')
1686!---------------------
1687END SUBROUTINE histend
1688!===
1689SUBROUTINE histwrite_r1d (pfileid,pvarname,pitau,pdata,nbindex,nindex)
1690!---------------------------------------------------------------------
1691  IMPLICIT NONE
1692!-
1693  INTEGER,INTENT(IN) :: pfileid, pitau, nbindex, nindex(nbindex)
1694  REAL,DIMENSION(:),INTENT(IN) :: pdata
1695  CHARACTER(LEN=*),INTENT(IN) :: pvarname
1696!-
1697  LOGICAL :: do_oper, do_write, largebuf
1698  INTEGER :: varid, io, nbpt_in, nbpt_out
1699  REAL, ALLOCATABLE,SAVE :: buff_tmp(:)
1700  INTEGER,SAVE :: buff_tmp_sz
1701  CHARACTER(LEN=7) :: tmp_opp
1702!-
1703  LOGICAL :: check = .FALSE.
1704!---------------------------------------------------------------------
1705!-
1706! 1.0 Try to catch errors like specifying the wrong file ID.
1707!     Thanks Marine for showing us what errors users can make !
1708!-
1709  IF ( (pfileid < 1).OR.(pfileid > nb_files) ) THEN
1710    CALL ipslerr (3,"histwrite", &
1711 &    'Illegal file ID in the histwrite of variable',pvarname,' ')
1712  ENDIF
1713!-
1714! 1.1 Find the id of the variable to be written and the real time
1715!-
1716  CALL histvar_seq (pfileid,pvarname,varid)
1717!-
1718! 2.0 do nothing for never operation
1719!-
1720  tmp_opp = topp(pfileid,varid)
1721!-
1722  IF (TRIM(tmp_opp) == "never") THEN
1723    last_opp_chk(pfileid,varid) = -99
1724    last_wrt_chk(pfileid,varid) = -99
1725  ENDIF
1726!-
1727! 3.0 We check if we need to do an operation
1728!-
1729  IF (last_opp_chk(pfileid,varid) == pitau) THEN
1730    CALL ipslerr (3,"histwrite", &
1731 &    'This variable as already been analysed at the present', &
1732 &    'time step',' ')
1733  ENDIF
1734!-
1735  CALL isittime &
1736 &  (pitau, date0(pfileid), deltat(pfileid), freq_opp(pfileid,varid), &
1737 &   last_opp(pfileid,varid), last_opp_chk(pfileid,varid), do_oper)
1738!-
1739! 4.0 We check if we need to write the data
1740!-
1741  IF (last_wrt_chk(pfileid,varid) == pitau) THEN
1742    CALL ipslerr (3,"histwrite", &
1743 &    'This variable as already been written for the present', &
1744 &    'time step',' ')
1745  ENDIF
1746!-
1747  CALL isittime &
1748 &  (pitau, date0(pfileid), deltat(pfileid), freq_wrt(pfileid,varid), &
1749 &   last_wrt(pfileid,varid), last_wrt_chk(pfileid,varid), do_write)
1750!-
1751! 5.0 histwrite called
1752!-
1753  IF (do_oper.OR.do_write) THEN
1754!-
1755!-- 5.1 Get the sizes of the data we will handle
1756!-
1757    IF (datasz_in(pfileid,varid,1) <= 0) THEN
1758!---- There is the risk here that the user has over-sized the array.
1759!---- But how can we catch this ?
1760!---- In the worst case we will do impossible operations
1761!---- on part of the data !
1762      datasz_in(pfileid,varid,1) = SIZE(pdata)
1763      datasz_in(pfileid,varid,2) = -1
1764      datasz_in(pfileid,varid,3) = -1
1765    ENDIF
1766!-
1767!-- 5.2 The maximum size of the data will give the size of the buffer
1768!-
1769    IF (datasz_max(pfileid,varid) <= 0) THEN
1770      largebuf = .FALSE.
1771      DO io=1,nbopp(pfileid,varid)
1772        IF (INDEX(fuchnbout,sopps(pfileid,varid,io)) > 0) THEN
1773          largebuf = .TRUE.
1774        ENDIF
1775      ENDDO
1776      IF (largebuf) THEN
1777        datasz_max(pfileid,varid) = &
1778 &        scsize(pfileid,varid,1) &
1779 &       *scsize(pfileid,varid,2) &
1780 &       *scsize(pfileid,varid,3)
1781      ELSE
1782        datasz_max(pfileid,varid) = &
1783 &        datasz_in(pfileid,varid,1)
1784      ENDIF
1785    ENDIF
1786!-
1787    IF (.NOT.ALLOCATED(buff_tmp)) THEN
1788      IF (check) THEN
1789        WRITE(*,*) &
1790 &       "histwrite_r1d : allocate buff_tmp for buff_sz = ", &
1791 &       datasz_max(pfileid,varid)
1792      ENDIF
1793      ALLOCATE (buff_tmp(datasz_max(pfileid,varid)))
1794      buff_tmp_sz = datasz_max(pfileid,varid)
1795    ELSE IF (datasz_max(pfileid,varid) > buff_tmp_sz) THEN
1796      IF (check) THEN
1797        WRITE(*,*) &
1798 &       "histwrite_r1d : re-allocate buff_tmp for buff_sz = ", &
1799 &       datasz_max(pfileid,varid)
1800      ENDIF
1801      DEALLOCATE (buff_tmp)
1802      ALLOCATE (buff_tmp(datasz_max(pfileid,varid)))
1803      buff_tmp_sz = datasz_max(pfileid,varid)
1804    ENDIF
1805!-
1806!-- We have to do the first operation anyway.
1807!-- Thus we do it here and change the ranke
1808!-- of the data at the same time. This should speed up things.
1809!-
1810    nbpt_in = datasz_in(pfileid,varid,1)
1811    nbpt_out = datasz_max(pfileid,varid)
1812    CALL mathop (sopps(pfileid,varid,1), nbpt_in, pdata, &
1813 &               missing_val, nbindex, nindex, &
1814 &               scal(pfileid,varid,1), nbpt_out, buff_tmp)
1815    CALL histwrite_real (pfileid, varid, pitau, nbpt_out, &
1816 &            buff_tmp, nbindex, nindex, do_oper, do_write)
1817  ENDIF
1818!-
1819! 6.0 Manage time steps
1820!-
1821  IF ((TRIM(tmp_opp) /= "once").AND.(TRIM(tmp_opp) /= "never")) THEN
1822    last_opp_chk(pfileid,varid) = pitau
1823    last_wrt_chk(pfileid,varid) = pitau
1824  ELSE
1825    last_opp_chk(pfileid,varid) = -99
1826    last_wrt_chk(pfileid,varid) = -99
1827  ENDIF
1828!---------------------------
1829END SUBROUTINE histwrite_r1d
1830!===
1831SUBROUTINE histwrite_r2d (pfileid,pvarname,pitau,pdata,nbindex,nindex)
1832!---------------------------------------------------------------------
1833  IMPLICIT NONE
1834!-
1835  INTEGER,INTENT(IN) :: pfileid, pitau, nbindex, nindex(nbindex)
1836  REAL,DIMENSION(:,:),INTENT(IN) :: pdata
1837  CHARACTER(LEN=*),INTENT(IN) :: pvarname
1838!-
1839  LOGICAL :: do_oper, do_write, largebuf
1840  INTEGER :: varid, io, nbpt_in(1:2), nbpt_out
1841  REAL, ALLOCATABLE,SAVE :: buff_tmp(:)
1842  INTEGER,SAVE :: buff_tmp_sz
1843  CHARACTER(LEN=7) :: tmp_opp
1844!-
1845  LOGICAL :: check = .FALSE.
1846!---------------------------------------------------------------------
1847!-
1848! 1.0 Try to catch errors like specifying the wrong file ID.
1849!     Thanks Marine for showing us what errors users can make !
1850!-
1851  IF ( (pfileid < 1).OR.(pfileid > nb_files) ) THEN
1852    CALL ipslerr (3,"histwrite", &
1853 &    'Illegal file ID in the histwrite of variable',pvarname,' ')
1854  ENDIF
1855!-
1856! 1.1 Find the id of the variable to be written and the real time
1857!-
1858  CALL histvar_seq (pfileid,pvarname,varid)
1859!-
1860! 2.0 do nothing for never operation
1861!-
1862  tmp_opp = topp(pfileid,varid)
1863!-
1864  IF (TRIM(tmp_opp) == "never") THEN
1865    last_opp_chk(pfileid,varid) = -99
1866    last_wrt_chk(pfileid,varid) = -99
1867  ENDIF
1868!-
1869! 3.0 We check if we need to do an operation
1870!-
1871  IF (last_opp_chk(pfileid,varid) == pitau) THEN
1872    CALL ipslerr (3,"histwrite", &
1873 &    'This variable as already been analysed at the present', &
1874 &    'time step',' ')
1875  ENDIF
1876!-
1877  CALL isittime &
1878 &  (pitau, date0(pfileid), deltat(pfileid), freq_opp(pfileid,varid), &
1879 &   last_opp(pfileid,varid), last_opp_chk(pfileid,varid), do_oper)
1880!-
1881! 4.0 We check if we need to write the data
1882!-
1883  IF (last_wrt_chk(pfileid,varid) == pitau) THEN
1884    CALL ipslerr (3,"histwrite", &
1885 &    'This variable as already been written for the present', &
1886 &    'time step',' ')
1887  ENDIF
1888!-
1889  CALL isittime &
1890 &  (pitau, date0(pfileid), deltat(pfileid), freq_wrt(pfileid,varid), &
1891 &   last_wrt(pfileid,varid), last_wrt_chk(pfileid,varid), do_write)
1892!-
1893! 5.0 histwrite called
1894!-
1895  IF (do_oper.OR.do_write) THEN
1896!-
1897!-- 5.1 Get the sizes of the data we will handle
1898!-
1899    IF (datasz_in(pfileid,varid,1) <= 0) THEN
1900!---- There is the risk here that the user has over-sized the array.
1901!---- But how can we catch this ?
1902!---- In the worst case we will do impossible operations
1903!---- on part of the data !
1904      datasz_in(pfileid,varid,1) = SIZE(pdata, DIM=1)
1905      datasz_in(pfileid,varid,2) = SIZE(pdata, DIM=2)
1906      datasz_in(pfileid,varid,3) = -1
1907    ENDIF
1908!-
1909!-- 5.2 The maximum size of the data will give the size of the buffer
1910!-
1911    IF (datasz_max(pfileid,varid) <= 0) THEN
1912      largebuf = .FALSE.
1913      DO io=1,nbopp(pfileid,varid)
1914        IF (INDEX(fuchnbout,sopps(pfileid,varid,io)) > 0) THEN
1915          largebuf = .TRUE.
1916        ENDIF
1917      ENDDO
1918      IF (largebuf) THEN
1919        datasz_max(pfileid,varid) = &
1920 &        scsize(pfileid,varid,1) &
1921 &       *scsize(pfileid,varid,2) &
1922 &       *scsize(pfileid,varid,3)
1923      ELSE
1924        datasz_max(pfileid,varid) = &
1925 &        datasz_in(pfileid,varid,1) &
1926 &       *datasz_in(pfileid,varid,2)
1927      ENDIF
1928    ENDIF
1929!-
1930    IF (.NOT.ALLOCATED(buff_tmp)) THEN
1931      IF (check) THEN
1932        WRITE(*,*) &
1933 &       "histwrite_r2d : allocate buff_tmp for buff_sz = ", &
1934 &       datasz_max(pfileid,varid)
1935      ENDIF
1936      ALLOCATE (buff_tmp(datasz_max(pfileid,varid)))
1937      buff_tmp_sz = datasz_max(pfileid,varid)
1938    ELSE IF (datasz_max(pfileid,varid) > buff_tmp_sz) THEN
1939      IF (check) THEN
1940        WRITE(*,*) &
1941 &       "histwrite_r2d : re-allocate buff_tmp for buff_sz = ", &
1942 &       datasz_max(pfileid,varid)
1943      ENDIF
1944      DEALLOCATE (buff_tmp)
1945      ALLOCATE (buff_tmp(datasz_max(pfileid,varid)))
1946      buff_tmp_sz = datasz_max(pfileid,varid)
1947    ENDIF
1948!-
1949!-- We have to do the first operation anyway.
1950!-- Thus we do it here and change the ranke
1951!-- of the data at the same time. This should speed up things.
1952!-
1953    nbpt_in(1:2) = datasz_in(pfileid,varid,1:2)
1954    nbpt_out = datasz_max(pfileid,varid)
1955    CALL mathop (sopps(pfileid,varid,1), nbpt_in, pdata, &
1956 &               missing_val, nbindex, nindex, &
1957 &               scal(pfileid,varid,1), nbpt_out, buff_tmp)
1958    CALL histwrite_real (pfileid, varid, pitau, nbpt_out, &
1959 &            buff_tmp, nbindex, nindex, do_oper, do_write)
1960  ENDIF
1961!-
1962! 6.0 Manage time steps
1963!-
1964  IF ((TRIM(tmp_opp) /= "once").AND.(TRIM(tmp_opp) /= "never")) THEN
1965    last_opp_chk(pfileid,varid) = pitau
1966    last_wrt_chk(pfileid,varid) = pitau
1967  ELSE
1968    last_opp_chk(pfileid,varid) = -99
1969    last_wrt_chk(pfileid,varid) = -99
1970  ENDIF
1971!---------------------------
1972END SUBROUTINE histwrite_r2d
1973!===
1974SUBROUTINE histwrite_r3d (pfileid,pvarname,pitau,pdata,nbindex,nindex)
1975!---------------------------------------------------------------------
1976  IMPLICIT NONE
1977!-
1978  INTEGER,INTENT(IN) :: pfileid, pitau, nbindex, nindex(nbindex)
1979  REAL,DIMENSION(:,:,:),INTENT(IN) :: pdata
1980  CHARACTER(LEN=*),INTENT(IN) :: pvarname
1981!-
1982  LOGICAL :: do_oper, do_write, largebuf
1983  INTEGER :: varid, io, nbpt_in(1:3), nbpt_out
1984  REAL, ALLOCATABLE,SAVE :: buff_tmp(:)
1985  INTEGER,SAVE :: buff_tmp_sz
1986  CHARACTER(LEN=7) :: tmp_opp
1987!-
1988  LOGICAL :: check = .FALSE.
1989!---------------------------------------------------------------------
1990!-
1991! 1.0 Try to catch errors like specifying the wrong file ID.
1992!     Thanks Marine for showing us what errors users can make !
1993!-
1994  IF ( (pfileid < 1).OR.(pfileid > nb_files) ) THEN
1995    CALL ipslerr (3,"histwrite", &
1996 &    'Illegal file ID in the histwrite of variable',pvarname,' ')
1997  ENDIF
1998!-
1999! 1.1 Find the id of the variable to be written and the real time
2000!-
2001  CALL histvar_seq (pfileid,pvarname,varid)
2002!-
2003! 2.0 do nothing for never operation
2004!-
2005  tmp_opp = topp(pfileid,varid)
2006!-
2007  IF (TRIM(tmp_opp) == "never") THEN
2008    last_opp_chk(pfileid,varid) = -99
2009    last_wrt_chk(pfileid,varid) = -99
2010  ENDIF
2011!-
2012! 3.0 We check if we need to do an operation
2013!-
2014  IF (last_opp_chk(pfileid,varid) == pitau) THEN
2015    CALL ipslerr (3,"histwrite", &
2016 &    'This variable as already been analysed at the present', &
2017 &    'time step',' ')
2018  ENDIF
2019!-
2020  CALL isittime &
2021 &  (pitau, date0(pfileid), deltat(pfileid), freq_opp(pfileid,varid), &
2022 &   last_opp(pfileid,varid), last_opp_chk(pfileid,varid), do_oper)
2023!-
2024! 4.0 We check if we need to write the data
2025!-
2026  IF (last_wrt_chk(pfileid,varid) == pitau) THEN
2027    CALL ipslerr (3,"histwrite", &
2028 &    'This variable as already been written for the present', &
2029 &    'time step',' ')
2030  ENDIF
2031!-
2032  CALL isittime &
2033 &  (pitau, date0(pfileid), deltat(pfileid), freq_wrt(pfileid,varid), &
2034 &   last_wrt(pfileid,varid), last_wrt_chk(pfileid,varid), do_write)
2035!-
2036! 5.0 histwrite called
2037!-
2038  IF (do_oper.OR.do_write) THEN
2039!-
2040!-- 5.1 Get the sizes of the data we will handle
2041!-
2042    IF (datasz_in(pfileid,varid,1) <= 0) THEN
2043!---- There is the risk here that the user has over-sized the array.
2044!---- But how can we catch this ?
2045!---- In the worst case we will do impossible operations
2046!---- on part of the data !
2047      datasz_in(pfileid,varid,1) = SIZE(pdata, DIM=1)
2048      datasz_in(pfileid,varid,2) = SIZE(pdata, DIM=2)
2049      datasz_in(pfileid,varid,3) = SIZE(pdata, DIM=3)
2050    ENDIF
2051!-
2052!-- 5.2 The maximum size of the data will give the size of the buffer
2053!-
2054    IF (datasz_max(pfileid,varid) <= 0) THEN
2055      largebuf = .FALSE.
2056      DO io =1,nbopp(pfileid,varid)
2057        IF (INDEX(fuchnbout,sopps(pfileid,varid,io)) > 0) THEN
2058          largebuf = .TRUE.
2059        ENDIF
2060      ENDDO
2061      IF (largebuf) THEN
2062        datasz_max(pfileid,varid) = &
2063 &        scsize(pfileid,varid,1) &
2064 &       *scsize(pfileid,varid,2) &
2065 &       *scsize(pfileid,varid,3)
2066      ELSE
2067        datasz_max(pfileid,varid) = &
2068 &        datasz_in(pfileid,varid,1) &
2069 &       *datasz_in(pfileid,varid,2) &
2070 &       *datasz_in(pfileid,varid,3)
2071      ENDIF
2072    ENDIF
2073!-
2074    IF (.NOT.ALLOCATED(buff_tmp)) THEN
2075      IF (check) THEN
2076        WRITE(*,*) &
2077 &       "histwrite_r1d : allocate buff_tmp for buff_sz = ", &
2078 &       datasz_max(pfileid,varid)
2079      ENDIF
2080      ALLOCATE (buff_tmp(datasz_max(pfileid,varid)))
2081      buff_tmp_sz = datasz_max(pfileid,varid)
2082    ELSE IF (datasz_max(pfileid,varid) > buff_tmp_sz) THEN
2083      IF (check) THEN
2084        WRITE(*,*) &
2085 &       "histwrite_r1d : re-allocate buff_tmp for buff_sz = ", &
2086 &       datasz_max(pfileid,varid)
2087      ENDIF
2088      DEALLOCATE (buff_tmp)
2089      ALLOCATE (buff_tmp(datasz_max(pfileid,varid)))
2090      buff_tmp_sz = datasz_max(pfileid,varid)
2091    ENDIF
2092!-
2093!-- We have to do the first operation anyway.
2094!-- Thus we do it here and change the ranke
2095!-- of the data at the same time. This should speed up things.
2096!-
2097    nbpt_in(1:3) = datasz_in(pfileid,varid,1:3)
2098    nbpt_out = datasz_max(pfileid,varid)
2099    CALL mathop (sopps(pfileid,varid,1), nbpt_in, pdata, &
2100 &               missing_val, nbindex, nindex, &
2101 &               scal(pfileid,varid,1), nbpt_out, buff_tmp)
2102    CALL histwrite_real (pfileid, varid, pitau, nbpt_out, &
2103 &            buff_tmp, nbindex, nindex, do_oper, do_write)
2104  ENDIF
2105!-
2106! 6.0 Manage time steps
2107!-
2108  IF ((TRIM(tmp_opp) /= "once").AND.(TRIM(tmp_opp) /= "never")) THEN
2109    last_opp_chk(pfileid,varid) = pitau
2110    last_wrt_chk(pfileid,varid) = pitau
2111  ELSE
2112    last_opp_chk(pfileid,varid) = -99
2113    last_wrt_chk(pfileid,varid) = -99
2114  ENDIF
2115!---------------------------
2116END SUBROUTINE histwrite_r3d
2117!===
2118SUBROUTINE histwrite_real &
2119 & (pfileid,varid,pitau,nbdpt,buff_tmp,nbindex,nindex,do_oper,do_write)
2120!---------------------------------------------------------------------
2121!- This subroutine is internal and does the calculations and writing
2122!- if needed. At a later stage it should be split into an operation
2123!- and writing subroutines.
2124!---------------------------------------------------------------------
2125  IMPLICIT NONE
2126!-
2127  INTEGER,INTENT(IN) :: pfileid,pitau,varid, &
2128 &                      nbindex,nindex(nbindex),nbdpt
2129  REAL,DIMENSION(:)  :: buff_tmp
2130  LOGICAL,INTENT(IN) :: do_oper,do_write
2131!-
2132  INTEGER :: tsz, ncid, ncvarid
2133  INTEGER :: i, iret, ipt, itax
2134  INTEGER :: io, nbin, nbout
2135  INTEGER,DIMENSION(4) :: corner, edges
2136  INTEGER :: itime
2137!-
2138  REAL :: rtime
2139  CHARACTER(LEN=7) :: tmp_opp
2140  REAL,ALLOCATABLE,SAVE :: buff_tmp2(:)
2141  INTEGER,SAVE          :: buff_tmp2_sz
2142  REAL,ALLOCATABLE,SAVE :: buffer_used(:)
2143  INTEGER,SAVE          :: buffer_sz
2144!-
2145  LOGICAL :: check = .FALSE.
2146!---------------------------------------------------------------------
2147  IF (check) THEN
2148    WRITE(*,*) "histwrite 0.0 :  VAR : ", name(pfileid,varid)
2149    WRITE(*,*) "histwrite 0.0 : nbindex, nindex :", &
2150    &  nbindex,nindex(1:MIN(3,nbindex)),'...',nindex(MAX(1,nbindex-3):nbindex)
2151  ENDIF
2152!-
2153! The sizes which can be encoutered
2154!-
2155  tsz = zsize(pfileid,varid,1)*zsize(pfileid,varid,2)*zsize(pfileid,varid,3)
2156!-
2157! 1.0 We allocate the memory needed to store the data between write
2158!     and the temporary space needed for operations.
2159!     We have to keep precedent buffer if needed
2160!-
2161  IF (.NOT. ALLOCATED(buffer)) THEN
2162    IF (check) WRITE(*,*) "histwrite_real 1.0 allocate buffer ",buff_pos
2163    ALLOCATE(buffer(buff_pos))
2164    buffer_sz = buff_pos
2165    buffer(:)=0.0
2166  ELSE IF (buffer_sz < buff_pos) THEN
2167    IF (check) WRITE(*,*) "histwrite_real 1.0.1 re-allocate buffer for ",buff_pos," instead of ",SIZE(buffer)
2168    IF (SUM(buffer)/=0.0) THEN
2169      IF (check) WRITE (*,*) ' histwrite : buffer has been used. We have to save it before re-allocating it '
2170      ALLOCATE (buffer_used(buffer_sz))
2171      buffer_used(:)=buffer(:)
2172      DEALLOCATE (buffer)
2173      ALLOCATE (buffer(buff_pos))
2174      buffer_sz = buff_pos
2175      buffer(:SIZE(buffer_used))=buffer_used
2176      DEALLOCATE (buffer_used)
2177    ELSE
2178      IF (check) WRITE (*,*) ' histwrite : buffer has not been used. We have just to re-allocating it '
2179      DEALLOCATE (buffer)
2180      ALLOCATE (buffer(buff_pos))
2181      buffer_sz = buff_pos
2182      buffer(:)=0.0
2183    ENDIF
2184  ENDIF
2185!-
2186! The buffers are only deallocated when more space is needed. This
2187! reduces the umber of allocates but increases memory needs.
2188!-
2189  IF (.NOT.ALLOCATED(buff_tmp2)) THEN
2190    IF (check) THEN
2191      WRITE(*,*) "histwrite_real 1.1 allocate buff_tmp2 ",SIZE(buff_tmp)
2192    ENDIF
2193    ALLOCATE (buff_tmp2(datasz_max(pfileid,varid)))
2194    buff_tmp2_sz = datasz_max(pfileid,varid)
2195  ELSE IF (datasz_max(pfileid,varid) > buff_tmp2_sz) THEN
2196    IF (check) THEN
2197      WRITE(*,*) "histwrite_real 1.2 re-allocate buff_tmp2 : ", &
2198     & SIZE(buff_tmp)," instead of ",SIZE(buff_tmp2)
2199    ENDIF
2200    DEALLOCATE (buff_tmp2)
2201    ALLOCATE (buff_tmp2(datasz_max(pfileid,varid)))
2202    buff_tmp2_sz = datasz_max(pfileid,varid)
2203  ENDIF
2204!-
2205  rtime = pitau * deltat(pfileid)
2206  tmp_opp = topp(pfileid,varid)
2207!-
2208! 3.0 Do the operations or transfer the slab of data into buff_tmp
2209!-
2210  IF (check) WRITE(*,*) "histwrite: 3.0", pfileid
2211!-
2212! 3.1 DO the Operations only if needed
2213!-
2214  IF (do_oper) THEN
2215    i = pfileid
2216    nbout = nbdpt
2217!-
2218!-- 3.4 We continue the sequence of operations
2219!--     we started in the interface routine
2220!-
2221    DO io = 2, nbopp(i,varid),2
2222      nbin = nbout
2223      nbout = datasz_max(i,varid)
2224      CALL mathop(sopps(i,varid,io),nbin,buff_tmp,missing_val, &
2225 &      nbindex,nindex,scal(i,varid,io),nbout,buff_tmp2)
2226      IF (check) THEN
2227        WRITE(*,*) &
2228 &       "histwrite: 3.4a nbout : ",nbin,nbout,sopps(i,varid,io)
2229      ENDIF
2230!-
2231      nbin = nbout
2232      nbout = datasz_max(i,varid)
2233      CALL mathop(sopps(i,varid,io+1),nbin,buff_tmp2,missing_val, &
2234 &      nbindex,nindex,scal(i,varid,io+1),nbout,buff_tmp)
2235      IF (check) THEN
2236        WRITE(*,*) &
2237 &       "histwrite: 3.4b nbout : ",nbin,nbout,sopps(i,varid,io+1)
2238      ENDIF
2239    ENDDO
2240!-
2241!   3.5 Zoom into the data
2242!-
2243    IF (check) THEN
2244      WRITE(*,*) &
2245 &     "histwrite: 3.5 size(buff_tmp) : ",SIZE(buff_tmp)
2246      WRITE(*,*) &
2247 &     "histwrite: 3.5 slab in X :",zorig(i,varid,1),zsize(i,varid,1)
2248      WRITE(*,*) &
2249 &     "histwrite: 3.5 slab in Y :",zorig(i,varid,2),zsize(i,varid,2)
2250      WRITE(*,*) &
2251 &     "histwrite: 3.5 slab in Z :",zorig(i,varid,3),zsize(i,varid,3)
2252      WRITE(*,*) &
2253 &     "histwrite: 3.5 slab of input:", &
2254 &     scsize(i,varid,1),scsize(i,varid,2),scsize(i,varid,3)
2255    ENDIF
2256    CALL trans_buff &
2257 &      (zorig(i,varid,1),zsize(i,varid,1), &
2258 &       zorig(i,varid,2),zsize(i,varid,2), &
2259 &       zorig(i,varid,3),zsize(i,varid,3), &
2260 &       scsize(i,varid,1),scsize(i,varid,2),scsize(i,varid,3), &
2261 &       buff_tmp, buff_tmp2_sz,buff_tmp2)
2262!-
2263!-- 4.0 Get the min and max of the field (buff_tmp)
2264!-
2265    IF (check) WRITE(*,*) "histwrite: 4.0 buff_tmp",pfileid,varid, &
2266 &    TRIM(tmp_opp),'----',LEN_TRIM(tmp_opp),nbindex
2267!-
2268    IF (hist_calc_rng(pfileid,varid)) THEN
2269      hist_minmax(pfileid,varid,1) = &
2270 &      MIN(hist_minmax(pfileid,varid,1), &
2271 &      MINVAL(buff_tmp2(1:tsz),MASK=buff_tmp2(1:tsz) /= missing_val))
2272      hist_minmax(pfileid,varid,2) = &
2273 &      MAX(hist_minmax(pfileid,varid,2), &
2274 &      MAXVAL(buff_tmp2(1:tsz),MASK=buff_tmp2(1:tsz) /= missing_val))
2275    ENDIF
2276!-
2277!-- 5.0 Do the operations if needed. In the case of instantaneous
2278!--     output we do not transfer to the buffer.
2279!-
2280    IF (check) WRITE(*,*) "histwrite: 5.0", pfileid, "tsz :", tsz
2281!-
2282    ipt = point(pfileid,varid)
2283!-
2284!   WRITE(*,*) 'OPE ipt, buffer :', pvarname, ipt, varid
2285!-
2286    IF (     (TRIM(tmp_opp) /= "inst") &
2287   &    .AND.(TRIM(tmp_opp) /= "once") ) THEN
2288      CALL moycum(tmp_opp,tsz,buffer(ipt:), &
2289     &       buff_tmp2,nb_opp(pfileid,varid))
2290    ENDIF
2291!-
2292    last_opp(pfileid,varid) = pitau
2293    nb_opp(pfileid,varid) = nb_opp(pfileid,varid)+1
2294!-
2295   ENDIF
2296!-
2297! 6.0 Write to file if needed
2298!-
2299  IF (check) WRITE(*,*) "histwrite: 6.0", pfileid
2300!-
2301  IF (do_write) THEN
2302!-
2303    ncvarid = ncvar_ids(pfileid,varid)
2304    ncid = ncdf_ids(pfileid)
2305!-
2306!-- 6.1 Do the operations that are needed before writting
2307!-
2308    IF (check) WRITE(*,*) "histwrite: 6.1", pfileid
2309!-
2310    IF (     (TRIM(tmp_opp) /= "inst") &
2311   &    .AND.(TRIM(tmp_opp) /= "once") ) THEN
2312      rtime = (rtime+last_wrt(pfileid,varid)*deltat(pfileid))/2.0
2313    ENDIF
2314!-
2315!-- 6.2 Add a value to the time axis of this variable if needed
2316!-
2317    IF (     (TRIM(tmp_opp) /= "l_max") &
2318   &    .AND.(TRIM(tmp_opp) /= "l_min") &
2319   &    .AND.(TRIM(tmp_opp) /= "once") ) THEN
2320!-
2321      IF (check) WRITE(*,*) "histwrite: 6.2", pfileid
2322!-
2323      itax = var_axid(pfileid,varid)
2324      itime = nb_wrt(pfileid,varid)+1
2325!-
2326      IF (tax_last(pfileid, itax) < itime) THEN
2327        iret = NF90_PUT_VAR (ncid,tdimid(pfileid,itax),(/ rtime /), &
2328 &                            start=(/ itime /),count=(/ 1 /))
2329        tax_last(pfileid, itax) = itime
2330      ENDIF
2331    ELSE
2332      itime=1
2333    ENDIF
2334!-
2335!-- 6.3 Write the data. Only in the case of instantaneous output
2336!       we do not write the buffer.
2337!-
2338    IF (check) THEN
2339      WRITE(*,*) "histwrite: 6.3",pfileid,ncid,ncvarid,varid,itime
2340    ENDIF
2341!-
2342    IF (scsize(pfileid,varid,3) == 1) THEN
2343      IF (regular(pfileid)) THEN
2344        corner(1:4) = (/ 1, 1, itime, 0 /)
2345        edges(1:4) = (/ zsize(pfileid,varid,1), &
2346 &                      zsize(pfileid,varid,2), &
2347 &                       1, 0 /)
2348      ELSE
2349        corner(1:4) = (/ 1, itime, 0, 0 /)
2350        edges(1:4) = (/ zsize(pfileid,varid,1), 1, 0, 0 /)
2351      ENDIF
2352    ELSE
2353      IF (regular(pfileid)) THEN
2354        corner(1:4) = (/ 1, 1, 1, itime /)
2355        edges(1:4) = (/ zsize(pfileid,varid,1), &
2356 &                      zsize(pfileid,varid,2), &
2357 &                      zsize(pfileid,varid,3), 1 /)
2358      ELSE
2359        corner(1:4) = (/ 1, 1, itime, 0 /)
2360        edges(1:4) = (/ zsize(pfileid,varid,1), &
2361 &                      zsize(pfileid,varid,3), 1, 0 /)
2362      ENDIF
2363    ENDIF
2364!-
2365    ipt = point(pfileid,varid)
2366!-
2367    IF (     (TRIM(tmp_opp) /= "inst") &
2368 &      .AND.(TRIM(tmp_opp) /= "once") ) THEN
2369      iret = NF90_PUT_VAR (ncid,ncvarid,buffer(ipt:), &
2370 &                       start=corner(1:4),count=edges(1:4))
2371    ELSE
2372      iret = NF90_PUT_VAR (ncid,ncvarid,buff_tmp2, &
2373 &                       start=corner(1:4),count=edges(1:4))
2374    ENDIF
2375!-
2376    last_wrt(pfileid,varid) = pitau
2377    nb_wrt(pfileid,varid) = nb_wrt(pfileid,varid)+1
2378    nb_opp(pfileid,varid) = 0
2379!---
2380!   After the write the file can be synchronized so that no data is
2381!   lost in case of a crash. This feature gives up on the benefits of
2382!   buffering and should only be used in debuging mode. A flag is
2383!   needed here to switch to this mode.
2384!---
2385!   iret = NF90_SYNC (ncid)
2386!-
2387  ENDIF
2388!----------------------------
2389END SUBROUTINE histwrite_real
2390!===
2391SUBROUTINE histvar_seq (pfid,pvarname,pvid)
2392!---------------------------------------------------------------------
2393!- This subroutine optimized the search for the variable in the table.
2394!- In a first phase it will learn the succession of the variables
2395!- called and then it will use the table to guess what comes next.
2396!- It is the best solution to avoid lengthy searches through array
2397!- vectors.
2398!-
2399!- ARGUMENTS :
2400!-
2401!- pfid  : id of the file on which we work
2402!- pvarname : The name of the variable we are looking for
2403!- pvid     : The var id we found
2404!---------------------------------------------------------------------
2405  IMPLICIT NONE
2406!-
2407  INTEGER,INTENT(in)  :: pfid
2408  CHARACTER(LEN=*),INTENT(IN) :: pvarname
2409  INTEGER,INTENT(out) :: pvid
2410!-
2411  LOGICAL,SAVE :: learning(nb_files_max)=.TRUE.
2412  INTEGER,SAVE :: overlap(nb_files_max) = -1
2413  INTEGER,SAVE :: varseq(nb_files_max,nb_var_max*3)
2414  INTEGER,SAVE :: varseq_len(nb_files_max) = 0
2415  INTEGER,SAVE :: varseq_pos(nb_files_max)
2416  INTEGER,SAVE :: varseq_err(nb_files_max) = 0
2417  INTEGER      :: ib,sp,nx,pos
2418  CHARACTER(LEN=70) :: str70
2419!-
2420  LOGICAL :: check = .FALSE.
2421!---------------------------------------------------------------------
2422  IF (check) THEN
2423    WRITE(*,*) 'histvar_seq, start of the subroutine :',learning(pfid)
2424  ENDIF
2425!-
2426  IF (learning(pfid)) THEN
2427!-
2428!-- 1.0 We compute the length over which we are going
2429!--     to check the overlap
2430!-
2431    IF (overlap(pfid) <= 0) THEN
2432      IF (nb_var(pfid) > 6) THEN
2433        overlap(pfid) = nb_var(pfid)/3*2
2434      ELSE
2435        overlap(pfid) = nb_var(pfid)
2436      ENDIF
2437    ENDIF
2438!-
2439!-- 1.1 Find the position of this string
2440!-
2441    CALL find_str (name(pfid,1:nb_var(pfid)),pvarname,pos)
2442    IF (pos > 0) THEN
2443      pvid = pos
2444    ELSE
2445      CALL ipslerr (3,"histvar_seq", &
2446 &      'The name of the variable you gave has not been declared', &
2447 &      'You should use subroutine histdef for declaring variable', &
2448 &      TRIM(pvarname))
2449    ENDIF
2450!-
2451!-- 1.2 If we have not given up we store the position
2452!--     in the sequence of calls
2453!-
2454    IF (varseq_err(pfid) >= 0) THEN
2455      sp = varseq_len(pfid)+1
2456      IF (sp <= nb_var_max*3) THEN
2457        varseq(pfid,sp) = pvid
2458        varseq_len(pfid) = sp
2459      ELSE
2460        CALL ipslerr (2,"histvar_seq",&
2461 &       'The learning process has failed and we give up. '// &
2462 &       'Either you sequence is',&
2463 &       'too complex or I am too dumb. '// &
2464 &       'This will only affect the efficiency',&
2465 &       'of your code. Thus if you wish to save time'// &
2466 &       ' contact the IOIPSL team. ')
2467        WRITE(*,*) 'The sequence we have found up to now :'
2468        WRITE(*,*) varseq(pfid,1:sp-1)
2469        varseq_err(pfid) = -1
2470      ENDIF
2471!-
2472!---- 1.3 Check if we have found the right overlap
2473!-
2474      IF (varseq_len(pfid) .GE. overlap(pfid)*2) THEN
2475!-
2476!------ We skip a few variables if needed as they could come
2477!------ from the initialisation of the model.
2478!-
2479        DO ib = 0, sp-overlap(pfid)*2
2480          IF ( learning(pfid) .AND.&
2481            & SUM(ABS(varseq(pfid,ib+1:ib+overlap(pfid)) -&
2482            & varseq(pfid,sp-overlap(pfid)+1:sp))) == 0 ) THEN
2483            learning(pfid) = .FALSE.
2484            varseq_len(pfid) = sp-overlap(pfid)-ib
2485            varseq_pos(pfid) = overlap(pfid)+ib
2486            varseq(pfid,1:varseq_len(pfid)) = &
2487 &            varseq(pfid,ib+1:ib+varseq_len(pfid))
2488          ENDIF
2489        ENDDO
2490      ENDIF
2491    ENDIF
2492  ELSE
2493!-
2494!-- 2.0 Now we know how the calls to histwrite are sequenced
2495!--     and we can get a guess at the var ID
2496!-
2497    nx = varseq_pos(pfid)+1
2498    IF (nx > varseq_len(pfid)) nx = 1
2499!-
2500    pvid = varseq(pfid,nx)
2501!-
2502    IF (TRIM(name(pfid,pvid)) /= TRIM(pvarname)) THEN
2503      CALL find_str (name(pfid,1:nb_var(pfid)),pvarname,pos)
2504      IF (pos > 0) THEN
2505        pvid = pos
2506      ELSE
2507        CALL ipslerr (3,"histvar_seq", &
2508 &  'The name of the variable you gave has not been declared',&
2509 &  'You should use subroutine histdef for declaring variable', &
2510 &  TRIM(pvarname))
2511      ENDIF
2512      varseq_err(pfid) = varseq_err(pfid)+1
2513    ELSE
2514!-
2515!---- We only keep the new position if we have found the variable
2516!---- this way. This way an out of sequence call to histwrite does
2517!---- not defeat the process.
2518!-
2519      varseq_pos(pfid) = nx
2520    ENDIF
2521!-
2522    IF (varseq_err(pfid) .GE. 10) THEN
2523      WRITE(str70,'("for file ",I3)') pfid
2524      CALL ipslerr (2,"histvar_seq", &
2525 &  'There were 10 errors in the learned sequence of variables',&
2526 &  str70,'This looks like a bug, please report it.')
2527         varseq_err(pfid) = 0
2528    ENDIF
2529  ENDIF
2530!-
2531  IF (check) THEN
2532    WRITE(*,*) &
2533 &   'histvar_seq, end of the subroutine :',TRIM(pvarname),pvid
2534  ENDIF
2535!-------------------------
2536END SUBROUTINE histvar_seq
2537!===
2538SUBROUTINE histsync (file)
2539!---------------------------------------------------------------------
2540!- This subroutine will synchronise all
2541!- (or one if defined) opened files.
2542!-
2543!- VERSION
2544!-
2545!---------------------------------------------------------------------
2546  IMPLICIT NONE
2547!-
2548! file  : optional argument for fileid
2549  INTEGER,INTENT(in),OPTIONAL :: file
2550!-
2551  INTEGER :: ifile,ncid,iret
2552!-
2553  LOGICAL :: file_exists
2554  LOGICAL :: check = .FALSE.
2555!---------------------------------------------------------------------
2556  IF (check) WRITE(*,*) 'Entering loop on files :', nb_files
2557!-
2558! 1.The loop on files to synchronise
2559!-
2560  DO ifile = 1,nb_files
2561!-
2562    IF (PRESENT(file)) THEN
2563      file_exists = (ifile == file)
2564    ELSE
2565      file_exists = .TRUE.
2566    ENDIF
2567!-
2568    IF (file_exists) THEN
2569      IF (check) THEN
2570        WRITE(*,*) 'Synchronising specified file number :', file
2571      ENDIF
2572      ncid = ncdf_ids(ifile)
2573      iret = NF90_SYNC (ncid)
2574    ENDIF
2575!-
2576  ENDDO
2577!----------------------
2578END SUBROUTINE histsync
2579!===
2580SUBROUTINE histclo (fid)
2581!---------------------------------------------------------------------
2582!- This subroutine will close all (or one if defined) opened files
2583!-
2584!- VERSION
2585!-
2586!---------------------------------------------------------------------
2587  IMPLICIT NONE
2588!-
2589! fid  : optional argument for fileid
2590  INTEGER,INTENT(in),OPTIONAL :: fid
2591!-
2592  INTEGER :: ifile,ncid,iret,iv
2593  INTEGER :: start_loop,end_loop
2594  CHARACTER(LEN=70) :: str70
2595!-
2596  LOGICAL :: check=.FALSE.
2597!---------------------------------------------------------------------
2598  IF (check) WRITE(*,*) 'Entering loop on files :', nb_files
2599!-
2600  IF (PRESENT(fid)) THEN
2601    start_loop = fid
2602    end_loop = fid
2603  ELSE
2604    start_loop = 1
2605    end_loop = nb_files
2606  ENDIF
2607!-
2608  DO ifile=start_loop,end_loop
2609    IF (check) WRITE(*,*) 'Closing specified file number :', ifile
2610    ncid = ncdf_ids(ifile)
2611    iret = NF90_REDEF (ncid)
2612!---
2613!-- 1. Loop on the number of variables to add some final information
2614!---
2615    IF (check) WRITE(*,*) 'Entering loop on vars :', nb_var(ifile)
2616    DO iv=1,nb_var(ifile)
2617      IF (hist_wrt_rng(ifile,iv)) THEN
2618        IF (check) THEN
2619          WRITE(*,*) 'min value for file :',ifile,' var n. :',iv, &
2620         &           ' is : ',hist_minmax(ifile,iv,1)
2621          WRITE(*,*) 'max value for file :',ifile,' var n. :',iv, &
2622         &           ' is : ',hist_minmax(ifile,iv,2)
2623        ENDIF
2624        IF (hist_calc_rng(ifile,iv)) THEN
2625!-------- Put the min and max values on the file
2626          iret = NF90_PUT_ATT (ncid,ncvar_ids(ifile,iv),'valid_min', &
2627 &                             REAL(hist_minmax(ifile,iv,1),KIND=4))
2628          iret = NF90_PUT_ATT (ncid,ncvar_ids(ifile,iv),'valid_max', &
2629 &                             REAL(hist_minmax(ifile,iv,2),KIND=4))
2630        ENDIF
2631      ENDIF
2632    ENDDO
2633!---
2634!-- 2. Close the file
2635!---
2636    IF (check) WRITE(*,*) 'close file :',ncid
2637    iret = NF90_CLOSE (ncid)
2638    IF (iret /= NF90_NOERR) THEN
2639      WRITE(str70,'("This file has been already closed :",I3)') ifile
2640      CALL ipslerr (2,'histclo',str70,'','')
2641    ENDIF
2642  ENDDO
2643!---------------------
2644END SUBROUTINE histclo
2645!===
2646SUBROUTINE ioconf_modname (str)
2647!---------------------------------------------------------------------
2648!- This subroutine allows to configure the name
2649!- of the model written into the file
2650!---------------------------------------------------------------------
2651  IMPLICIT NONE
2652!-
2653  CHARACTER(LEN=*),INTENT(IN) :: str
2654!---------------------------------------------------------------------
2655  IF (.NOT.lock_modname) THEN
2656    model_name = str(1:MIN(LEN_TRIM(str),80))
2657    lock_modname = .TRUE.
2658  ELSE
2659    CALL ipslerr (2,"ioconf_modname", &
2660   &  'The model name can only be changed once and only', &
2661   &  'before it is used. It is now set to :',model_name)
2662  ENDIF
2663!----------------------------
2664END SUBROUTINE ioconf_modname
2665!-
2666!===
2667!-
2668!-----------------
2669END MODULE histcom
Note: See TracBrowser for help on using the repository browser.