source: branches/UKMO/r6232_tracer_advection/NEMOGCM/EXTERNAL/IOIPSL/src/histcom.f90 @ 9295

Last change on this file since 9295 was 9295, checked in by jcastill, 3 years ago

Remove svn keywords

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