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

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