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

Last change on this file since 964 was 964, checked in by bellier, 15 years ago

Reinitialize some elements at the closure of a file

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