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

Last change on this file since 1023 was 1023, checked in by bellier, 12 years ago

New version with bounds for reductive time operations

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