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

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

The positive direction of a vertical axis must be "up" or "down".
If the optional argument pdirect of histvert is not given or is invalid,
the attribute will not be written in the Netcdf file.

  • 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=20) :: direction
826  INTEGER :: iret,leng,nfid
827  LOGICAL :: l_dbg
828!---------------------------------------------------------------------
829  CALL ipsldbg (old_status=l_dbg)
830!-
831! 1.0 Verifications :
832!    Do we have enough space for an extra axis ?
833!    Is the name already in use ?
834!-
835  IF (l_dbg) WRITE(*,*) "histvert : 1.0 Verifications", &
836 &                      pzaxname,'---',pzaxunit,'---',pzaxtitle
837!-
838! Direction of the vertical axis. Can we get if from the user.
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 (     PRESENT(pdirect)    &
850 &    .AND.(direction /= 'up') &
851 &    .AND.(direction /= 'down') ) THEN
852    direction = 'unknown'
853    CALL ipslerr (2,"histvert",&
854 & "The specified positive direction for the vertical axis is invalid.",&
855 & "The value must be up or down.","The attribute will not be written.")
856  ENDIF
857!-
858  IF (W_F(idf)%n_zax+1 > nb_zax_max) THEN
859    CALL ipslerr (3,"histvert", &
860   &  'Table of vertical axes too small. You should increase ',&
861   &  'nb_zax_max in histcom.f90 in order to accomodate all ', &
862   &  'these variables ')
863  ENDIF
864!-
865  iv = W_F(idf)%n_zax
866  IF (iv > 1) THEN
867    CALL find_str (W_F(idf)%zax_name(1:iv-1),pzaxname,pos)
868  ELSE
869    pos = 0
870  ENDIF
871!-
872  IF (pos > 0) THEN
873    WRITE(str71,'("Check variable ",A," in file",I3)') &
874 &    TRIM(pzaxname),idf
875    CALL ipslerr (3,"histvert", &
876 &    "Vertical axis already exists",TRIM(str71), &
877 &    "Can also be a wrong file ID in another declaration")
878  ENDIF
879!-
880  iv = W_F(idf)%n_zax+1
881!-
882! 2.0 Add the information to the file
883!-
884  IF (l_dbg) &
885 &  WRITE(*,*) "histvert : 2.0 Add the information to the file"
886!-
887  nfid = W_F(idf)%ncfid
888!-
889  leng = MIN(LEN_TRIM(pzaxname),20)
890  iret = NF90_DEF_DIM (nfid,pzaxname(1:leng),pzsize,zaxid_tmp)
891  iret = NF90_DEF_VAR (nfid,pzaxname(1:leng),NF90_REAL4, &
892 &                     zaxid_tmp,zdimid)
893  iret = NF90_PUT_ATT (nfid,zdimid,'axis',"Z")
894  iret = NF90_PUT_ATT (nfid,zdimid,'standard_name',"model_level_number")
895  leng = MIN(LEN_TRIM(pzaxunit),20)
896  IF (leng > 0) THEN
897    iret = NF90_PUT_ATT (nfid,zdimid,'units',pzaxunit(1:leng))
898  ENDIF
899  IF (direction /= 'unknown') THEN
900    iret = NF90_PUT_ATT (nfid,zdimid,'positive',TRIM(direction))
901  ENDIF
902  iret = NF90_PUT_ATT (nfid,zdimid,'valid_min', &
903 &                     REAL(MINVAL(pzvalues(1:pzsize)),KIND=4))
904  iret = NF90_PUT_ATT (nfid,zdimid,'valid_max', &
905 &                     REAL(MAXVAL(pzvalues(1:pzsize)),KIND=4))
906  leng = MIN(LEN_TRIM(pzaxname),20)
907  iret = NF90_PUT_ATT (nfid,zdimid,'title',pzaxname(1:leng))
908  leng = MIN(LEN_TRIM(pzaxtitle),80)
909  iret = NF90_PUT_ATT (nfid,zdimid,'long_name',pzaxtitle(1:leng))
910!-
911  iret = NF90_ENDDEF (nfid)
912!-
913  iret = NF90_PUT_VAR (nfid,zdimid,pzvalues(1:pzsize))
914!-
915  iret = NF90_REDEF (nfid)
916!-
917!- 3.0 add the information to the common
918!-
919  IF (l_dbg) &
920  &  WRITE(*,*) "histvert : 3.0 add the information to the common"
921!-
922  W_F(idf)%n_zax = iv
923  W_F(idf)%zax_size(iv) = pzsize
924  W_F(idf)%zax_name(iv) = pzaxname
925  W_F(idf)%zax_ids(iv) = zaxid_tmp
926  pzaxid = iv
927!----------------------
928END SUBROUTINE histvert
929!===
930SUBROUTINE histdef &
931 &  (idf,pvarname,ptitle,punit, &
932 &   pxsize,pysize,phoriid,pzsize,par_oriz,par_szz,pzid, &
933 &   xtype,popp,pfreq_opp,pfreq_wrt,var_range,standard_name)
934!---------------------------------------------------------------------
935!- With this subroutine each variable to be archived on the history
936!- tape should be declared.
937!-
938!- It gives the user the choise of operation
939!- to be performed on the variables, the frequency of this operation
940!- and finaly the frequency of the archiving.
941!-
942!- INPUT
943!-
944!- idf      : ID of the file the variable should be archived in
945!- pvarname : Name of the variable, short and easy to remember
946!- ptitle   : Full name of the variable
947!- punit    : Units of the variable (no units if blank string)
948!-
949!- The next 3 arguments give the size of that data
950!- that will be passed to histwrite. The zoom will be
951!- done there with the horizontal information obtained
952!- in histbeg and the vertical information to follow.
953!-
954!- pxsize   : Size in X direction (size of the data that will be
955!-            given to histwrite)
956!- pysize   : Size in Y direction
957!- phoriid  : ID of the horizontal axis
958!-
959!- The next two arguments give the vertical zoom to use.
960!-
961!- pzsize   : Size in Z direction (If 1 then no axis is declared
962!-            for this variable and pzid is not used)
963!- par_oriz : Off set of the zoom
964!- par_szz  : Size of the zoom
965!-
966!- pzid     : ID of the vertical axis to use. It has to have
967!-            the size of the zoom.
968!- xtype    : External netCDF type (hist_r4/hist_r8)
969!- popp     : Operation to be performed. The following options
970!-            exist today :
971!- inst : keeps instantaneous values for writting
972!- ave  : Computes the average from call between writes
973!- pfreq_opp: Frequency of this operation (in seconds)
974!- pfreq_wrt: Frequency at which the variable should be
975!-            written (in seconds)
976!- var_range: Range of the variable.
977!-            If the minimum is greater than the maximum,
978!-            the values will be calculated.
979!-
980!- VERSION
981!---------------------------------------------------------------------
982  IMPLICIT NONE
983!-
984  INTEGER,INTENT(IN) :: idf,pxsize,pysize,pzsize,pzid
985  INTEGER,INTENT(IN) :: par_oriz,par_szz,xtype,phoriid
986  CHARACTER(LEN=*),INTENT(IN) :: pvarname,punit,popp,ptitle
987  REAL,INTENT(IN) :: pfreq_opp,pfreq_wrt
988  REAL,DIMENSION(2),OPTIONAL,INTENT(IN) :: var_range
989  CHARACTER(LEN=*),OPTIONAL,INTENT(IN) :: standard_name
990!-
991  INTEGER :: iv
992  CHARACTER(LEN=70) :: str70,str71,str72
993  CHARACTER(LEN=20) :: tmp_name
994  CHARACTER(LEN=40) :: str40
995  CHARACTER(LEN=10) :: str10
996  CHARACTER(LEN=120) :: ex_topps
997  REAL :: un_an,un_jour,test_fopp,test_fwrt
998  INTEGER :: pos,buff_sz
999  LOGICAL :: l_dbg
1000!---------------------------------------------------------------------
1001  CALL ipsldbg (old_status=l_dbg)
1002!-
1003  ex_topps = 'ave, inst, t_min, t_max, t_sum, once, never, l_max, l_min'
1004!-
1005  W_F(idf)%n_var = W_F(idf)%n_var+1
1006  iv = W_F(idf)%n_var
1007!-
1008  IF (iv > nb_var_max) THEN
1009    CALL ipslerr (3,"histdef", &
1010   &  'Table of variables too small. You should increase nb_var_max',&
1011   &  'in histcom.f90 in order to accomodate all these variables', &
1012   &  ' ')
1013  ENDIF
1014!-
1015! 1.0 Transfer informations on the variable to the common
1016!     and verify that it does not already exist
1017!-
1018  IF (l_dbg) WRITE(*,*) "histdef : 1.0"
1019!-
1020  IF (iv > 1) THEN
1021    CALL find_str (W_F(idf)%W_V(1:iv-1)%v_name,pvarname,pos)
1022  ELSE
1023    pos = 0
1024  ENDIF
1025!-
1026  IF (pos > 0) THEN
1027    str70 = "Variable already exists"
1028    WRITE(str71,'("Check variable  ",a," in file",I3)') &
1029 &    TRIM(pvarname),idf
1030    str72 = "Can also be a wrong file ID in another declaration"
1031    CALL ipslerr (3,"histdef",str70,str71,str72)
1032  ENDIF
1033!-
1034  W_F(idf)%W_V(iv)%v_name = pvarname
1035  W_F(idf)%W_V(iv)%title = ptitle
1036  W_F(idf)%W_V(iv)%unit_name = punit
1037  IF (PRESENT(standard_name)) THEN
1038    W_F(idf)%W_V(iv)%std_name = standard_name
1039  ELSE
1040    W_F(idf)%W_V(iv)%std_name = ptitle
1041  ENDIF
1042  tmp_name = W_F(idf)%W_V(iv)%v_name
1043!-
1044! 1.1 decode the operations
1045!-
1046  W_F(idf)%W_V(iv)%fullop = popp
1047  CALL buildop &
1048 &  (TRIM(popp),ex_topps,W_F(idf)%W_V(iv)%topp,missing_val, &
1049 &   W_F(idf)%W_V(iv)%sopp,W_F(idf)%W_V(iv)%scal, &
1050 &   W_F(idf)%W_V(iv)%nbopp)
1051!-
1052! 1.2 If we have an even number of operations
1053!     then we need to add identity
1054!-
1055  IF ( MOD(W_F(idf)%W_V(iv)%nbopp,2) == 0) THEN
1056    W_F(idf)%W_V(iv)%nbopp = W_F(idf)%W_V(iv)%nbopp+1
1057    W_F(idf)%W_V(iv)%sopp(W_F(idf)%W_V(iv)%nbopp) = 'ident'
1058    W_F(idf)%W_V(iv)%scal(W_F(idf)%W_V(iv)%nbopp) = missing_val
1059  ENDIF
1060!-
1061! 1.3 External type of the variable
1062!-
1063  IF (xtype == hist_r8) THEN
1064    W_F(idf)%W_V(iv)%v_typ = hist_r8
1065  ELSE
1066    W_F(idf)%W_V(iv)%v_typ = hist_r4
1067  ENDIF
1068!-
1069! 2.0 Put the size of the variable in the common and check
1070!-
1071  IF (l_dbg) THEN
1072    WRITE(*,*) "histdef : 2.0",idf,iv,W_F(idf)%W_V(iv)%nbopp, &
1073 &    W_F(idf)%W_V(iv)%sopp(1:W_F(idf)%W_V(iv)%nbopp), &
1074 &    W_F(idf)%W_V(iv)%scal(1:W_F(idf)%W_V(iv)%nbopp)
1075  ENDIF
1076!-
1077  W_F(idf)%W_V(iv)%scsize(1:3) = (/ pxsize,pysize,pzsize /)
1078  W_F(idf)%W_V(iv)%zorig(1:3) = &
1079 &  (/ W_F(idf)%slab_ori(1),W_F(idf)%slab_ori(2),par_oriz /)
1080  W_F(idf)%W_V(iv)%zsize(1:3) = &
1081 &  (/ W_F(idf)%slab_siz(1),W_F(idf)%slab_siz(2),par_szz /)
1082!-
1083! Is the size of the full array the same as that of the coordinates ?
1084!-
1085  IF (    (pxsize > W_F(idf)%full_size(1)) &
1086 &    .OR.(pysize > W_F(idf)%full_size(2)) ) THEN
1087!-
1088    str70 = "The size of the variable is different "// &
1089 &          "from the one of the coordinates"
1090    WRITE(str71,'("Size of coordinates :",2I4)') &
1091 &   W_F(idf)%full_size(1),W_F(idf)%full_size(2)
1092    WRITE(str72,'("Size declared for variable ",a," :",2I4)') &
1093 &   TRIM(tmp_name),pxsize,pysize
1094    CALL ipslerr (3,"histdef",str70,str71,str72)
1095  ENDIF
1096!-
1097! Is the size of the zoom smaller than the coordinates ?
1098!-
1099  IF (    (W_F(idf)%full_size(1) < W_F(idf)%slab_siz(1)) &
1100 &    .OR.(W_F(idf)%full_size(2) < W_F(idf)%slab_siz(2)) ) THEN
1101    str70 = &
1102 &   "Size of variable should be greater or equal to those of the zoom"
1103    WRITE(str71,'("Size of XY zoom :",2I4)') &
1104 &   W_F(idf)%slab_siz(1),W_F(idf)%slab_siz(2)
1105    WRITE(str72,'("Size declared for variable ",A," :",2I4)') &
1106 &   TRIM(tmp_name),pxsize,pysize
1107    CALL ipslerr (3,"histdef",str70,str71,str72)
1108  ENDIF
1109!-
1110! 2.1 We store the horizontal grid information with minimal
1111!     and a fall back onto the default grid
1112!-
1113  IF ( (phoriid > 0).AND.(phoriid <= W_F(idf)%n_hax) ) THEN
1114    W_F(idf)%W_V(iv)%h_axid = phoriid
1115  ELSE
1116    W_F(idf)%W_V(iv)%h_axid = 1
1117    CALL ipslerr (2,"histdef", &
1118   &  'We use the default grid for variable as an invalide',&
1119   &  'ID was provided for variable : ',TRIM(pvarname))
1120  ENDIF
1121!-
1122! 2.2 Check the vertical coordinates if needed
1123!-
1124  IF (par_szz > 1) THEN
1125!-
1126!-- Does the vertical coordinate exist ?
1127!-
1128    IF (pzid > W_F(idf)%n_zax) THEN
1129      WRITE(str70, &
1130 &    '("The vertical coordinate chosen for variable ",A)') &
1131 &     TRIM(tmp_name)
1132      str71 = " Does not exist."
1133      CALL ipslerr (3,"histdef",str70,str71," ")
1134    ENDIF
1135!-
1136!-- Is the vertical size of the variable equal to that of the axis ?
1137!-
1138    IF (par_szz /= W_F(idf)%zax_size(pzid)) THEN
1139      str70 = "The size of the zoom does not correspond "// &
1140 &            "to the size of the chosen vertical axis"
1141      WRITE(str71,'("Size of zoom in z :",I4)') par_szz
1142      WRITE(str72,'("Size declared for axis ",A," :",I4)') &
1143 &     TRIM(W_F(idf)%zax_name(pzid)),W_F(idf)%zax_size(pzid)
1144      CALL ipslerr (3,"histdef",str70,str71,str72)
1145    ENDIF
1146!-
1147!-- Is the zoom smaller that the total size of the variable ?
1148!-
1149    IF (pzsize < par_szz) THEN
1150      str70 = "The vertical size of variable "// &
1151 &            "is smaller than that of the zoom."
1152      WRITE(str71,'("Declared vertical size of data :",I5)') pzsize
1153      WRITE(str72,'("Size of zoom for variable ",a," = ",I5)') &
1154 &     TRIM(tmp_name),par_szz
1155      CALL ipslerr (3,"histdef",str70,str71,str72)
1156    ENDIF
1157    W_F(idf)%W_V(iv)%z_axid = pzid
1158  ELSE
1159    W_F(idf)%W_V(iv)%z_axid = -99
1160  ENDIF
1161!-
1162! 3.0 We get the size of the arrays histwrite will get
1163!     and eventually allocate the time_buffer
1164!-
1165  IF (l_dbg) THEN
1166    WRITE(*,*) "histdef : 3.0"
1167  ENDIF
1168!-
1169  buff_sz = W_F(idf)%W_V(iv)%zsize(1) &
1170 &         *W_F(idf)%W_V(iv)%zsize(2) &
1171 &         *W_F(idf)%W_V(iv)%zsize(3)
1172!-
1173  IF (     (TRIM(W_F(idf)%W_V(iv)%topp) /= "inst") &
1174 &    .AND.(TRIM(W_F(idf)%W_V(iv)%topp) /= "once") &
1175 &    .AND.(TRIM(W_F(idf)%W_V(iv)%topp) /= "never") )THEN
1176    ALLOCATE(W_F(idf)%W_V(iv)%t_bf(buff_sz))
1177    W_F(idf)%W_V(iv)%t_bf(:) = 0.
1178    IF (l_dbg) THEN
1179      WRITE(*,*) "histdef : 3.0 allocating time_buffer for", &
1180 &      " idf = ",idf," iv = ",iv," size = ",buff_sz
1181    ENDIF
1182  ENDIF
1183!-
1184! 4.0 Transfer the frequency of the operations and check
1185!     for validity. We have to pay attention to negative values
1186!     of the frequency which indicate monthly time-steps.
1187!     The strategy is to bring it back to seconds for the tests
1188!-
1189  IF (l_dbg) WRITE(*,*) "histdef : 4.0"
1190!-
1191  W_F(idf)%W_V(iv)%freq_opp = pfreq_opp
1192  W_F(idf)%W_V(iv)%freq_wrt = pfreq_wrt
1193!-
1194  CALL ioget_calendar(un_an,un_jour)
1195  IF (pfreq_opp < 0) THEN
1196    CALL ioget_calendar(un_an)
1197    test_fopp = pfreq_opp*(-1.)*un_an/12.*un_jour
1198  ELSE
1199    test_fopp = pfreq_opp
1200  ENDIF
1201  IF (pfreq_wrt < 0) THEN
1202    CALL ioget_calendar(un_an)
1203    test_fwrt = pfreq_wrt*(-1.)*un_an/12.*un_jour
1204  ELSE
1205    test_fwrt = pfreq_wrt
1206  ENDIF
1207!-
1208! 4.1 Frequency of operations and output should be larger than deltat !
1209!-
1210  IF (test_fopp < W_F(idf)%deltat) THEN
1211    str70 = 'Frequency of operations should be larger than deltat'
1212    WRITE(str71,'("It is not the case for variable ",a," :",F10.4)') &
1213 &   TRIM(tmp_name),pfreq_opp
1214    str72 = "PATCH : frequency set to deltat"
1215!-
1216    CALL ipslerr (2,"histdef",str70,str71,str72)
1217!-
1218    W_F(idf)%W_V(iv)%freq_opp = W_F(idf)%deltat
1219  ENDIF
1220!-
1221  IF (test_fwrt < W_F(idf)%deltat) THEN
1222    str70 = 'Frequency of output should be larger than deltat'
1223    WRITE(str71,'("It is not the case for variable ",a," :",F10.4)') &
1224 &   TRIM(tmp_name),pfreq_wrt
1225    str72 = "PATCH : frequency set to deltat"
1226!-
1227    CALL ipslerr (2,"histdef",str70,str71,str72)
1228!-
1229    W_F(idf)%W_V(iv)%freq_wrt = W_F(idf)%deltat
1230  ENDIF
1231!-
1232! 4.2 First the existence of the operation is tested and then
1233!     its compaticility with the choice of frequencies
1234!-
1235  IF (TRIM(W_F(idf)%W_V(iv)%topp) == "inst") THEN
1236    IF (test_fopp /= test_fwrt) THEN
1237      str70 = 'For instantaneous output the frequency '// &
1238 &            'of operations and output'
1239      WRITE(str71, &
1240 &     '("should be the same, this was not case for variable ",a)') &
1241 &      TRIM(tmp_name)
1242      str72 = "PATCH : The smalest frequency of both is used"
1243      CALL ipslerr (2,"histdef",str70,str71,str72)
1244      IF (test_fopp < test_fwrt) THEN
1245        W_F(idf)%W_V(iv)%freq_opp = pfreq_opp
1246        W_F(idf)%W_V(iv)%freq_wrt = pfreq_opp
1247      ELSE
1248        W_F(idf)%W_V(iv)%freq_opp = pfreq_wrt
1249        W_F(idf)%W_V(iv)%freq_wrt = pfreq_wrt
1250      ENDIF
1251    ENDIF
1252  ELSE IF (INDEX(ex_topps,TRIM(W_F(idf)%W_V(iv)%topp)) > 0) THEN
1253    IF (test_fopp > test_fwrt) THEN
1254      str70 = 'For averages the frequency of operations '// &
1255 &            'should be smaller or equal'
1256      WRITE(str71, &
1257 &     '("to that of output. It is not the case for variable ",a)') &
1258 &     TRIM(tmp_name)
1259      str72 = 'PATCH : The output frequency is used for both'
1260      CALL ipslerr (2,"histdef",str70,str71,str72)
1261      W_F(idf)%W_V(iv)%freq_opp = pfreq_wrt
1262    ENDIF
1263  ELSE
1264    WRITE (str70,'("Operation on variable ",A," is unknown")') &
1265 &   TRIM(tmp_name)
1266    WRITE (str71,'("operation requested is :",A)') &
1267 &   W_F(idf)%W_V(iv)%topp
1268    WRITE (str72,'("File ID :",I3)') idf
1269    CALL ipslerr (3,"histdef",str70,str71,str72)
1270  ENDIF
1271!-
1272! 5.0 Initialize other variables of the common
1273!-
1274  IF (l_dbg) WRITE(*,*) "histdef : 5.0"
1275!-
1276  W_F(idf)%W_V(iv)%hist_wrt_rng = (PRESENT(var_range))
1277  IF (W_F(idf)%W_V(iv)%hist_wrt_rng) THEN
1278    W_F(idf)%W_V(iv)%hist_calc_rng = (var_range(1) > var_range(2))
1279    IF (W_F(idf)%W_V(iv)%hist_calc_rng) THEN
1280      W_F(idf)%W_V(iv)%hist_minmax(1:2) = &
1281 &      (/ ABS(missing_val),-ABS(missing_val) /)
1282    ELSE
1283      W_F(idf)%W_V(iv)%hist_minmax(1:2) = var_range(1:2)
1284    ENDIF
1285  ENDIF
1286!-
1287! - freq_opp(idf,iv)/2./deltat(idf)
1288  W_F(idf)%W_V(iv)%last_opp = W_F(idf)%itau0
1289! - freq_wrt(idf,iv)/2./deltat(idf)
1290  W_F(idf)%W_V(iv)%last_wrt = W_F(idf)%itau0
1291! - freq_opp(idf,iv)/2./deltat(idf)
1292  W_F(idf)%W_V(iv)%last_opp_chk = W_F(idf)%itau0
1293! - freq_wrt(idf,iv)/2./deltat(idf)
1294  W_F(idf)%W_V(iv)%last_wrt_chk = W_F(idf)%itau0
1295  W_F(idf)%W_V(iv)%nb_opp = 0
1296  W_F(idf)%W_V(iv)%nb_wrt = 0
1297!-
1298! 6.0 Get the time axis for this variable
1299!-
1300  IF (l_dbg) WRITE(*,*) "histdef : 6.0"
1301!-
1302! No time axis for once, l_max, l_min or never operation
1303!-
1304  IF (     (TRIM(W_F(idf)%W_V(iv)%topp) /= 'once')  &
1305 &    .AND.(TRIM(W_F(idf)%W_V(iv)%topp) /= 'never') &
1306 &    .AND.(TRIM(W_F(idf)%W_V(iv)%topp) /= 'l_max') &
1307 &    .AND.(TRIM(W_F(idf)%W_V(iv)%topp) /= 'l_min') ) THEN
1308    IF (TRIM(W_F(idf)%W_V(iv)%topp) == 'inst') THEN
1309      str10 = 't_inst_'
1310    ELSE
1311      str10 = 't_op_'
1312    ENDIF
1313    IF (W_F(idf)%W_V(iv)%freq_wrt > 0) THEN
1314      WRITE (UNIT=str40,FMT='(A,I8.8)') &
1315&      TRIM(str10),INT(W_F(idf)%W_V(iv)%freq_wrt)
1316    ELSE
1317      WRITE (UNIT=str40,FMT='(A,I2.2,"month")') &
1318&      TRIM(str10),ABS(INT(W_F(idf)%W_V(iv)%freq_wrt))
1319    ENDIF
1320    CALL find_str (W_F(idf)%W_V(1:W_F(idf)%n_tax)%tax_name,str40,pos)
1321    IF (pos < 0) THEN
1322      W_F(idf)%n_tax = W_F(idf)%n_tax+1
1323      W_F(idf)%W_V(iv)%l_bnd = &
1324 &      (TRIM(W_F(idf)%W_V(iv)%topp) /= 'inst')
1325      W_F(idf)%W_V(W_F(idf)%n_tax)%tax_name = str40
1326      W_F(idf)%W_V(W_F(idf)%n_tax)%tax_last = 0
1327      W_F(idf)%W_V(iv)%t_axid = W_F(idf)%n_tax
1328    ELSE
1329      W_F(idf)%W_V(iv)%t_axid = pos
1330    ENDIF
1331  ELSE
1332    IF (l_dbg) THEN
1333      WRITE(*,*) "histdef : 7.0 ",TRIM(W_F(idf)%W_V(iv)%topp),'----'
1334    ENDIF
1335    W_F(idf)%W_V(iv)%t_axid = -99
1336  ENDIF
1337!-
1338! 7.0 prepare frequence of writing and operation
1339!     for never or once operation
1340!-
1341  IF (    (TRIM(W_F(idf)%W_V(iv)%topp) == 'once')  &
1342 &    .OR.(TRIM(W_F(idf)%W_V(iv)%topp) == 'never') ) THEN
1343    W_F(idf)%W_V(iv)%freq_opp = 0.
1344    W_F(idf)%W_V(iv)%freq_wrt = 0.
1345  ENDIF
1346!---------------------
1347END SUBROUTINE histdef
1348!===
1349SUBROUTINE histend (idf)
1350!---------------------------------------------------------------------
1351!- This subroutine end the decalaration of variables and sets the
1352!- time axes in the netcdf file and puts it into the write mode.
1353!-
1354!- INPUT
1355!-
1356!- idf : ID of the file to be worked on
1357!-
1358!- VERSION
1359!-
1360!---------------------------------------------------------------------
1361  IMPLICIT NONE
1362!-
1363  INTEGER,INTENT(IN) :: idf
1364!-
1365  INTEGER :: nfid,nvid,iret,ndim,iv,itx,ziv,itax,dim_cnt
1366  INTEGER,DIMENSION(4) :: dims
1367  INTEGER :: year,month,day,hours,minutes
1368  REAL :: sec
1369  REAL :: rtime0
1370  CHARACTER(LEN=30) :: str30
1371  CHARACTER(LEN=35) :: str35
1372  CHARACTER(LEN=120) :: assoc
1373  CHARACTER(LEN=70) :: str70
1374  CHARACTER(LEN=3),DIMENSION(12) :: cal =   &
1375 &  (/ 'JAN','FEB','MAR','APR','MAY','JUN', &
1376 &     'JUL','AUG','SEP','OCT','NOV','DEC' /)
1377  CHARACTER(LEN=7) :: tmp_opp
1378  LOGICAL :: l_b
1379  LOGICAL :: l_dbg
1380!---------------------------------------------------------------------
1381  CALL ipsldbg (old_status=l_dbg)
1382!-
1383  nfid = W_F(idf)%ncfid
1384!-
1385! 1.0 Create the time axes
1386!-
1387  IF (l_dbg) WRITE(*,*) "histend : 1.0"
1388!-
1389! 1.1 Define the time dimensions needed for this file
1390!-
1391  iret = NF90_DEF_DIM (nfid,'time_counter', &
1392 &                     NF90_UNLIMITED,W_F(idf)%tid)
1393  DO iv=1,W_F(idf)%n_var
1394    IF (W_F(idf)%W_V(iv)%l_bnd) THEN
1395      iret = NF90_DEF_DIM (nfid,'tbnds',2,W_F(idf)%bid)
1396      EXIT
1397    ENDIF
1398  ENDDO
1399!-
1400! 1.2 Define all the time axes needed for this file
1401!-
1402  DO itx=1,W_F(idf)%n_tax
1403    dims(1) = W_F(idf)%tid
1404    l_b = (INDEX(W_F(idf)%W_V(itx)%tax_name,"t_op_") == 1)
1405    IF (itx > 1) THEN
1406      str30 = W_F(idf)%W_V(itx)%tax_name
1407    ELSE
1408      str30 = "time_counter"
1409    ENDIF
1410    IF (l_b) THEN
1411      str35 = TRIM(str30)//'_bnds'
1412    ENDIF
1413    iret = NF90_DEF_VAR (nfid,TRIM(str30),NF90_REAL8, &
1414 &          dims(1),W_F(idf)%W_V(itx)%tdimid)
1415    IF (itx <= 1) THEN
1416      iret = NF90_PUT_ATT (nfid,W_F(idf)%W_V(itx)%tdimid,'axis',"T")
1417    ENDIF
1418    iret = NF90_PUT_ATT (nfid,W_F(idf)%W_V(itx)%tdimid, &
1419 &          'standard_name',"time")
1420!---
1421!   To transform the current itau into a real date and take it
1422!   as the origin of the file requires the time counter to change.
1423!   Thus it is an operation the user has to ask for.
1424!   This function should thus only be re-instated
1425!   if there is a ioconf routine to control it.
1426!---
1427!-- rtime0 = itau2date(itau0(idf),date0(idf),deltat(idf))
1428    rtime0 = W_F(idf)%date0
1429!-
1430    CALL ju2ymds(rtime0,year,month,day,sec)
1431!---
1432!   Catch any error induced by a change in calendar !
1433!---
1434    IF (year < 0) THEN
1435      year = 2000+year
1436    ENDIF
1437!-
1438    hours = INT(sec/(60.*60.))
1439    minutes = INT((sec-hours*60.*60.)/60.)
1440    sec = sec-(hours*60.*60.+minutes*60.)
1441!-
1442    WRITE (UNIT=str70, &
1443 &   FMT='(A,I4.4,"-",I2.2,"-",I2.2," ",I2.2,":",I2.2,":",I2.2)') &
1444 &    'seconds since ',year,month,day,hours,minutes,INT(sec)
1445    iret = NF90_PUT_ATT (nfid,W_F(idf)%W_V(itx)%tdimid, &
1446 &           'units',TRIM(str70))
1447!-
1448    CALL ioget_calendar (str30)
1449    iret = NF90_PUT_ATT (nfid,W_F(idf)%W_V(itx)%tdimid, &
1450 &           'calendar',TRIM(str30))
1451!-
1452    iret = NF90_PUT_ATT (nfid,W_F(idf)%W_V(itx)%tdimid, &
1453 &           'title','Time')
1454!-
1455    iret = NF90_PUT_ATT (nfid,W_F(idf)%W_V(itx)%tdimid, &
1456 &           'long_name','Time axis')
1457!-
1458    WRITE (UNIT=str70, &
1459 &   FMT='(" ",I4.4,"-",A3,"-",I2.2," ",I2.2,":",I2.2,":",I2.2)') &
1460 &    year,cal(month),day,hours,minutes,INT(sec)
1461    iret = NF90_PUT_ATT (nfid,W_F(idf)%W_V(itx)%tdimid, &
1462 &           'time_origin',TRIM(str70))
1463!---
1464    IF (l_b) THEN
1465      iret = NF90_PUT_ATT (nfid,W_F(idf)%W_V(itx)%tdimid, &
1466 &             'bounds',TRIM(str35))
1467      dims(1:2) = (/ W_F(idf)%bid,W_F(idf)%tid /)
1468      iret = NF90_DEF_VAR (nfid,TRIM(str35),NF90_REAL8, &
1469 &            dims(1:2),W_F(idf)%W_V(itx)%tbndid)
1470    ENDIF
1471  ENDDO
1472!-
1473! 2.0 declare the variables
1474!-
1475  IF (l_dbg) WRITE(*,*) "histend : 2.0"
1476!-
1477  DO iv=1,W_F(idf)%n_var
1478!---
1479    itax = W_F(idf)%W_V(iv)%t_axid
1480!---
1481    IF (W_F(idf)%regular) THEN
1482      dims(1:2) = (/ W_F(idf)%xid,W_F(idf)%yid /)
1483      dim_cnt = 2
1484    ELSE
1485      dims(1) = W_F(idf)%xid
1486      dim_cnt = 1
1487    ENDIF
1488!---
1489    tmp_opp = W_F(idf)%W_V(iv)%topp
1490    ziv = W_F(idf)%W_V(iv)%z_axid
1491!---
1492!   2.1 dimension of field
1493!---
1494    IF ((TRIM(tmp_opp) /= 'never')) THEN
1495      IF (     (TRIM(tmp_opp) /= 'once')  &
1496     &    .AND.(TRIM(tmp_opp) /= 'l_max') &
1497     &    .AND.(TRIM(tmp_opp) /= 'l_min') ) THEN
1498        IF (ziv == -99) THEN
1499          ndim = dim_cnt+1
1500          dims(dim_cnt+1:dim_cnt+2) = (/ W_F(idf)%tid,0 /)
1501        ELSE
1502          ndim = dim_cnt+2
1503          dims(dim_cnt+1:dim_cnt+2) = &
1504 &         (/ W_F(idf)%zax_ids(ziv),W_F(idf)%tid /)
1505        ENDIF
1506      ELSE
1507        IF (ziv == -99) THEN
1508          ndim = dim_cnt
1509          dims(dim_cnt+1:dim_cnt+2) = (/ 0,0 /)
1510        ELSE
1511          ndim = dim_cnt+1
1512          dims(dim_cnt+1:dim_cnt+2) = (/ W_F(idf)%zax_ids(ziv),0 /)
1513        ENDIF
1514      ENDIF
1515!-
1516      iret = NF90_DEF_VAR (nfid,TRIM(W_F(idf)%W_V(iv)%v_name), &
1517 &             W_F(idf)%W_V(iv)%v_typ,dims(1:ABS(ndim)),nvid)
1518!-
1519      W_F(idf)%W_V(iv)%ncvid = nvid
1520!-
1521      IF (LEN_TRIM(W_F(idf)%W_V(iv)%unit_name) > 0) THEN
1522        iret = NF90_PUT_ATT (nfid,nvid,'units', &
1523 &                           TRIM(W_F(idf)%W_V(iv)%unit_name))
1524      ENDIF
1525      iret = NF90_PUT_ATT (nfid,nvid,'standard_name', &
1526 &                         TRIM(W_F(idf)%W_V(iv)%std_name))
1527!-
1528      IF (W_F(idf)%W_V(iv)%v_typ == hist_r8) THEN
1529        iret = NF90_PUT_ATT (nfid,nvid,'_FillValue',NF90_FILL_REAL8)
1530      ELSE
1531        iret = NF90_PUT_ATT (nfid,nvid,'_FillValue',NF90_FILL_REAL4)
1532      ENDIF
1533      IF (W_F(idf)%W_V(iv)%hist_wrt_rng) THEN
1534        IF (W_F(idf)%W_V(iv)%v_typ == hist_r8) THEN
1535          iret = NF90_PUT_ATT (nfid,nvid,'valid_min', &
1536 &                 REAL(W_F(idf)%W_V(iv)%hist_minmax(1),KIND=8))
1537          iret = NF90_PUT_ATT (nfid,nvid,'valid_max', &
1538 &                 REAL(W_F(idf)%W_V(iv)%hist_minmax(2),KIND=8))
1539        ELSE
1540          iret = NF90_PUT_ATT (nfid,nvid,'valid_min', &
1541 &                 REAL(W_F(idf)%W_V(iv)%hist_minmax(1),KIND=4))
1542          iret = NF90_PUT_ATT (nfid,nvid,'valid_max', &
1543 &                 REAL(W_F(idf)%W_V(iv)%hist_minmax(2),KIND=4))
1544        ENDIF
1545      ENDIF
1546      iret = NF90_PUT_ATT (nfid,nvid,'long_name', &
1547 &                         TRIM(W_F(idf)%W_V(iv)%title))
1548      iret = NF90_PUT_ATT (nfid,nvid,'online_operation', &
1549 &                         TRIM(W_F(idf)%W_V(iv)%fullop))
1550!-
1551      SELECT CASE(ndim)
1552      CASE(-3,2:4)
1553      CASE DEFAULT
1554        CALL ipslerr (3,"histend", &
1555       &  'less than 2 or more than 4 dimensions are not', &
1556       &  'allowed at this stage',' ')
1557      END SELECT
1558!-
1559      assoc=TRIM(W_F(idf)%hax_name(W_F(idf)%W_V(iv)%h_axid,2)) &
1560 &   //' '//TRIM(W_F(idf)%hax_name(W_F(idf)%W_V(iv)%h_axid,1))
1561!-
1562      ziv = W_F(idf)%W_V(iv)%z_axid
1563      IF (ziv > 0) THEN
1564        str30 = W_F(idf)%zax_name(ziv)
1565        assoc = TRIM(str30)//' '//TRIM(assoc)
1566      ENDIF
1567!-
1568      IF (itax > 0) THEN
1569        IF (itax > 1) THEN
1570          str30 = W_F(idf)%W_V(itax)%tax_name
1571        ELSE
1572          str30 = "time_counter"
1573        ENDIF
1574        assoc = TRIM(str30)//' '//TRIM(assoc)
1575!-
1576        IF (l_dbg) THEN
1577          WRITE(*,*) "histend : 2.0.n, freq_opp, freq_wrt", &
1578 &          W_F(idf)%W_V(iv)%freq_opp,W_F(idf)%W_V(iv)%freq_wrt
1579        ENDIF
1580!-
1581        iret = NF90_PUT_ATT (nfid,nvid,'interval_operation', &
1582 &                           REAL(W_F(idf)%W_V(iv)%freq_opp,KIND=4))
1583        iret = NF90_PUT_ATT (nfid,nvid,'interval_write', &
1584 &                           REAL(W_F(idf)%W_V(iv)%freq_wrt,KIND=4))
1585      ENDIF
1586      iret = NF90_PUT_ATT (nfid,nvid,'coordinates',TRIM(assoc))
1587    ENDIF
1588  ENDDO
1589!-
1590! 2.2 Add DOMAIN attributes if needed
1591!-
1592  IF (W_F(idf)%dom_id_svg >= 0) THEN
1593    CALL flio_dom_att (nfid,W_F(idf)%dom_id_svg)
1594  ENDIF
1595!-
1596! 3.0 Put the netcdf file into write mode
1597!-
1598  IF (l_dbg) WRITE(*,*) "histend : 3.0"
1599!-
1600  iret = NF90_ENDDEF (nfid)
1601!-
1602! 4.0 Give some informations to the user
1603!-
1604  IF (l_dbg) WRITE(*,*) "histend : 4.0"
1605!-
1606  WRITE(str70,'("All variables have been initialized on file :",I3)') idf
1607  CALL ipslerr (1,'histend',str70,'',' ')
1608!---------------------
1609END SUBROUTINE histend
1610!===
1611SUBROUTINE histwrite_r1d (idf,pvarname,pitau,pdata,nbindex,nindex)
1612!---------------------------------------------------------------------
1613  IMPLICIT NONE
1614!-
1615  INTEGER,INTENT(IN) :: idf,pitau,nbindex
1616  INTEGER,DIMENSION(nbindex),INTENT(IN) :: nindex
1617  REAL,DIMENSION(:),INTENT(IN) :: pdata
1618  CHARACTER(LEN=*),INTENT(IN) :: pvarname
1619!---------------------------------------------------------------------
1620  CALL histw_rnd (idf,pvarname,pitau,nbindex,nindex,pdata_1d=pdata)
1621!---------------------------
1622END SUBROUTINE histwrite_r1d
1623!===
1624SUBROUTINE histwrite_r2d (idf,pvarname,pitau,pdata,nbindex,nindex)
1625!---------------------------------------------------------------------
1626  IMPLICIT NONE
1627!-
1628  INTEGER,INTENT(IN) :: idf,pitau,nbindex
1629  INTEGER,DIMENSION(nbindex),INTENT(IN) :: nindex
1630  REAL,DIMENSION(:,:),INTENT(IN) :: pdata
1631  CHARACTER(LEN=*),INTENT(IN) :: pvarname
1632!---------------------------------------------------------------------
1633  CALL histw_rnd (idf,pvarname,pitau,nbindex,nindex,pdata_2d=pdata)
1634!---------------------------
1635END SUBROUTINE histwrite_r2d
1636!===
1637SUBROUTINE histwrite_r3d (idf,pvarname,pitau,pdata,nbindex,nindex)
1638!---------------------------------------------------------------------
1639  IMPLICIT NONE
1640!-
1641  INTEGER,INTENT(IN) :: idf,pitau,nbindex
1642  INTEGER,DIMENSION(nbindex),INTENT(IN) :: nindex
1643  REAL,DIMENSION(:,:,:),INTENT(IN) :: pdata
1644  CHARACTER(LEN=*),INTENT(IN) :: pvarname
1645!---------------------------------------------------------------------
1646  CALL histw_rnd (idf,pvarname,pitau,nbindex,nindex,pdata_3d=pdata)
1647!---------------------------
1648END SUBROUTINE histwrite_r3d
1649!===
1650SUBROUTINE histw_rnd (idf,pvarname,pitau,nbindex,nindex, &
1651  &                   pdata_1d,pdata_2d,pdata_3d)
1652!---------------------------------------------------------------------
1653  IMPLICIT NONE
1654!-
1655  INTEGER,INTENT(IN) :: idf,pitau,nbindex
1656  INTEGER,DIMENSION(nbindex),INTENT(IN) :: nindex
1657  CHARACTER(LEN=*),INTENT(IN) :: pvarname
1658  REAL,DIMENSION(:),INTENT(IN),OPTIONAL     :: pdata_1d
1659  REAL,DIMENSION(:,:),INTENT(IN),OPTIONAL   :: pdata_2d
1660  REAL,DIMENSION(:,:,:),INTENT(IN),OPTIONAL :: pdata_3d
1661!-
1662  LOGICAL :: do_oper,do_write,largebuf,l1d,l2d,l3d
1663  INTEGER :: iv,io,nbpt_out
1664  INTEGER              :: nbpt_in1
1665  INTEGER,DIMENSION(2) :: nbpt_in2
1666  INTEGER,DIMENSION(3) :: nbpt_in3
1667  REAL,ALLOCATABLE,DIMENSION(:),SAVE :: tbf_1
1668  CHARACTER(LEN=7) :: tmp_opp
1669  CHARACTER(LEN=13) :: c_nam
1670  LOGICAL :: l_dbg
1671!---------------------------------------------------------------------
1672  CALL ipsldbg (old_status=l_dbg)
1673!-
1674  l1d=PRESENT(pdata_1d); l2d=PRESENT(pdata_2d); l3d=PRESENT(pdata_3d);
1675  IF      (l1d) THEN
1676    c_nam = 'histwrite_r1d'
1677  ELSE IF (l2d) THEN
1678    c_nam = 'histwrite_r2d'
1679  ELSE IF (l3d) THEN
1680    c_nam = 'histwrite_r3d'
1681  ENDIF
1682!-
1683  IF (l_dbg) THEN
1684    WRITE(*,*) "histwrite : ",c_nam
1685  ENDIF
1686!-
1687! 1.0 Try to catch errors like specifying the wrong file ID.
1688!     Thanks Marine for showing us what errors users can make !
1689!-
1690  IF ( (idf < 1).OR.(idf > nb_files_max) ) THEN
1691    CALL ipslerr (3,"histwrite", &
1692 &    'Illegal file ID in the histwrite of variable',pvarname,' ')
1693  ENDIF
1694!-
1695! 1.1 Find the id of the variable to be written and the real time
1696!-
1697  CALL histvar_seq (idf,pvarname,iv)
1698!-
1699! 2.0 do nothing for never operation
1700!-
1701  tmp_opp = W_F(idf)%W_V(iv)%topp
1702!-
1703  IF (TRIM(tmp_opp) == "never") THEN
1704    W_F(idf)%W_V(iv)%last_opp_chk = -99
1705    W_F(idf)%W_V(iv)%last_wrt_chk = -99
1706  ENDIF
1707!-
1708! 3.0 We check if we need to do an operation
1709!-
1710  IF (W_F(idf)%W_V(iv)%last_opp_chk == pitau) THEN
1711    CALL ipslerr (3,"histwrite", &
1712 &    'This variable has already been analysed at the present', &
1713 &    'time step',TRIM(pvarname))
1714  ENDIF
1715!-
1716  CALL isittime &
1717 &  (pitau,W_F(idf)%date0,W_F(idf)%deltat, &
1718 &   W_F(idf)%W_V(iv)%freq_opp, &
1719 &   W_F(idf)%W_V(iv)%last_opp, &
1720 &   W_F(idf)%W_V(iv)%last_opp_chk,do_oper)
1721!-
1722! 4.0 We check if we need to write the data
1723!-
1724  IF (W_F(idf)%W_V(iv)%last_wrt_chk == pitau) THEN
1725    CALL ipslerr (3,"histwrite", &
1726 &    'This variable as already been written for the present', &
1727 &    'time step',' ')
1728  ENDIF
1729!-
1730  CALL isittime &
1731 &  (pitau,W_F(idf)%date0,W_F(idf)%deltat, &
1732 &   W_F(idf)%W_V(iv)%freq_wrt, &
1733 &   W_F(idf)%W_V(iv)%last_wrt, &
1734 &   W_F(idf)%W_V(iv)%last_wrt_chk,do_write)
1735!-
1736! 5.0 histwrite called
1737!-
1738  IF (do_oper.OR.do_write) THEN
1739!-
1740!-- 5.1 Get the sizes of the data we will handle
1741!-
1742    IF (W_F(idf)%W_V(iv)%datasz_in(1) <= 0) THEN
1743!---- There is the risk here that the user has over-sized the array.
1744!---- But how can we catch this ?
1745!---- In the worst case we will do impossible operations
1746!---- on part of the data !
1747      W_F(idf)%W_V(iv)%datasz_in(1:3) = -1
1748      IF      (l1d) THEN
1749        W_F(idf)%W_V(iv)%datasz_in(1) = SIZE(pdata_1d)
1750      ELSE IF (l2d) THEN
1751        W_F(idf)%W_V(iv)%datasz_in(1) = SIZE(pdata_2d,DIM=1)
1752        W_F(idf)%W_V(iv)%datasz_in(2) = SIZE(pdata_2d,DIM=2)
1753      ELSE IF (l3d) THEN
1754        W_F(idf)%W_V(iv)%datasz_in(1) = SIZE(pdata_3d,DIM=1)
1755        W_F(idf)%W_V(iv)%datasz_in(2) = SIZE(pdata_3d,DIM=2)
1756        W_F(idf)%W_V(iv)%datasz_in(3) = SIZE(pdata_3d,DIM=3)
1757      ENDIF
1758    ENDIF
1759!-
1760!-- 5.2 The maximum size of the data will give the size of the buffer
1761!-
1762    IF (W_F(idf)%W_V(iv)%datasz_max <= 0) THEN
1763      largebuf = .FALSE.
1764      DO io=1,W_F(idf)%W_V(iv)%nbopp
1765        IF (INDEX(fuchnbout,W_F(idf)%W_V(iv)%sopp(io)) > 0) THEN
1766          largebuf = .TRUE.
1767        ENDIF
1768      ENDDO
1769      IF (largebuf) THEN
1770        W_F(idf)%W_V(iv)%datasz_max = &
1771 &        W_F(idf)%W_V(iv)%scsize(1) &
1772 &       *W_F(idf)%W_V(iv)%scsize(2) &
1773 &       *W_F(idf)%W_V(iv)%scsize(3)
1774      ELSE
1775        IF      (l1d) THEN
1776          W_F(idf)%W_V(iv)%datasz_max = &
1777 &          W_F(idf)%W_V(iv)%datasz_in(1)
1778        ELSE IF (l2d) THEN
1779          W_F(idf)%W_V(iv)%datasz_max = &
1780 &          W_F(idf)%W_V(iv)%datasz_in(1) &
1781 &         *W_F(idf)%W_V(iv)%datasz_in(2)
1782        ELSE IF (l3d) THEN
1783          W_F(idf)%W_V(iv)%datasz_max = &
1784 &          W_F(idf)%W_V(iv)%datasz_in(1) &
1785 &         *W_F(idf)%W_V(iv)%datasz_in(2) &
1786 &         *W_F(idf)%W_V(iv)%datasz_in(3)
1787        ENDIF
1788      ENDIF
1789    ENDIF
1790!-
1791    IF (.NOT.ALLOCATED(tbf_1)) THEN
1792      IF (l_dbg) THEN
1793        WRITE(*,*) &
1794 &       c_nam//" : allocate tbf_1 for size = ", &
1795 &       W_F(idf)%W_V(iv)%datasz_max
1796      ENDIF
1797      ALLOCATE(tbf_1(W_F(idf)%W_V(iv)%datasz_max))
1798    ELSE IF (W_F(idf)%W_V(iv)%datasz_max > SIZE(tbf_1)) THEN
1799      IF (l_dbg) THEN
1800        WRITE(*,*) &
1801 &       c_nam//" : re-allocate tbf_1 for size = ", &
1802 &       W_F(idf)%W_V(iv)%datasz_max
1803      ENDIF
1804      DEALLOCATE(tbf_1)
1805      ALLOCATE(tbf_1(W_F(idf)%W_V(iv)%datasz_max))
1806    ENDIF
1807!-
1808!-- We have to do the first operation anyway.
1809!-- Thus we do it here and change the ranke
1810!-- of the data at the same time. This should speed up things.
1811!-
1812    nbpt_out = W_F(idf)%W_V(iv)%datasz_max
1813    IF      (l1d) THEN
1814      nbpt_in1 = W_F(idf)%W_V(iv)%datasz_in(1)
1815      CALL mathop (W_F(idf)%W_V(iv)%sopp(1),nbpt_in1,pdata_1d, &
1816 &                 missing_val,nbindex,nindex, &
1817 &                 W_F(idf)%W_V(iv)%scal(1),nbpt_out,tbf_1)
1818    ELSE IF (l2d) THEN
1819      nbpt_in2(1:2) = W_F(idf)%W_V(iv)%datasz_in(1:2)
1820      CALL mathop (W_F(idf)%W_V(iv)%sopp(1),nbpt_in2,pdata_2d, &
1821 &                 missing_val,nbindex,nindex, &
1822 &                 W_F(idf)%W_V(iv)%scal(1),nbpt_out,tbf_1)
1823    ELSE IF (l3d) THEN
1824      nbpt_in3(1:3) = W_F(idf)%W_V(iv)%datasz_in(1:3)
1825      CALL mathop (W_F(idf)%W_V(iv)%sopp(1),nbpt_in3,pdata_3d, &
1826 &                 missing_val,nbindex,nindex, &
1827 &                 W_F(idf)%W_V(iv)%scal(1),nbpt_out,tbf_1)
1828    ENDIF
1829    CALL histwrite_real (idf,iv,pitau,nbpt_out, &
1830 &            tbf_1,nbindex,nindex,do_oper,do_write)
1831  ENDIF
1832!-
1833! 6.0 Manage time steps
1834!-
1835  IF ((TRIM(tmp_opp) /= "once").AND.(TRIM(tmp_opp) /= "never")) THEN
1836    W_F(idf)%W_V(iv)%last_opp_chk = pitau
1837    W_F(idf)%W_V(iv)%last_wrt_chk = pitau
1838  ELSE
1839    W_F(idf)%W_V(iv)%last_opp_chk = -99
1840    W_F(idf)%W_V(iv)%last_wrt_chk = -99
1841  ENDIF
1842!-----------------------
1843END SUBROUTINE histw_rnd
1844!===
1845SUBROUTINE histwrite_real &
1846 & (idf,iv,pitau,nbdpt,tbf_1,nbindex,nindex,do_oper,do_write)
1847!---------------------------------------------------------------------
1848!- This subroutine is internal and does the calculations and writing
1849!- if needed. At a later stage it should be split into an operation
1850!- and writing subroutines.
1851!---------------------------------------------------------------------
1852  IMPLICIT NONE
1853!-
1854  INTEGER,INTENT(IN) :: idf,pitau,iv, &
1855 &                      nbindex,nindex(nbindex),nbdpt
1856  REAL,DIMENSION(:)  :: tbf_1
1857  LOGICAL,INTENT(IN) :: do_oper,do_write
1858!-
1859  INTEGER :: tsz,nfid,nvid,iret,itax,io,nbin,nbout
1860  INTEGER :: nx,ny,nz,ky,kz,kt,kc
1861  INTEGER,DIMENSION(4) :: corner,edges
1862  INTEGER :: itime
1863!-
1864  REAL :: rtime
1865  REAL,DIMENSION(2) :: t_bnd
1866  CHARACTER(LEN=7) :: tmp_opp
1867  REAL,ALLOCATABLE,DIMENSION(:),SAVE :: tbf_2
1868  LOGICAL :: l_dbg
1869!---------------------------------------------------------------------
1870  CALL ipsldbg (old_status=l_dbg)
1871!-
1872  IF (l_dbg) THEN
1873    WRITE(*,*) "histwrite 0.0 : VAR : ",W_F(idf)%W_V(iv)%v_name
1874    WRITE(*,*) "histwrite 0.0 : nbindex :",nbindex
1875    WRITE(*,*) "histwrite 0.0 : nindex  :",nindex(1:MIN(3,nbindex)),'...'
1876  ENDIF
1877!-
1878! The sizes which can be encoutered
1879!-
1880  tsz =  W_F(idf)%W_V(iv)%zsize(1) &
1881 &      *W_F(idf)%W_V(iv)%zsize(2) &
1882 &      *W_F(idf)%W_V(iv)%zsize(3)
1883!-
1884! 1.0 We allocate and the temporary space needed for operations.
1885!     The buffers are only deallocated when more space is needed.
1886!     This reduces the umber of allocates but increases memory needs.
1887!-
1888  IF (.NOT.ALLOCATED(tbf_2)) THEN
1889    IF (l_dbg) THEN
1890      WRITE(*,*) "histwrite_real 1.1 allocate tbf_2 ",SIZE(tbf_1)
1891    ENDIF
1892    ALLOCATE(tbf_2(W_F(idf)%W_V(iv)%datasz_max))
1893  ELSE IF (W_F(idf)%W_V(iv)%datasz_max > SIZE(tbf_2)) THEN
1894    IF (l_dbg) THEN
1895      WRITE(*,*) "histwrite_real 1.2 re-allocate tbf_2 : ", &
1896     & SIZE(tbf_1)," instead of ",SIZE(tbf_2)
1897    ENDIF
1898    DEALLOCATE(tbf_2)
1899    ALLOCATE(tbf_2(W_F(idf)%W_V(iv)%datasz_max))
1900  ENDIF
1901!-
1902  rtime = pitau*W_F(idf)%deltat
1903  tmp_opp = W_F(idf)%W_V(iv)%topp
1904!-
1905! 3.0 Do the operations or transfer the slab of data into tbf_1
1906!-
1907  IF (l_dbg) THEN
1908    WRITE(*,*) "histwrite: 3.0",idf
1909  ENDIF
1910!-
1911! 3.1 DO the Operations only if needed
1912!-
1913  IF (do_oper) THEN
1914    nbout = nbdpt
1915!-
1916!-- 3.4 We continue the sequence of operations
1917!--     we started in the interface routine
1918!-
1919    DO io=2,W_F(idf)%W_V(iv)%nbopp,2
1920      nbin = nbout
1921      nbout = W_F(idf)%W_V(iv)%datasz_max
1922      CALL mathop(W_F(idf)%W_V(iv)%sopp(io),nbin,tbf_1, &
1923 &      missing_val,nbindex,nindex,W_F(idf)%W_V(iv)%scal(io), &
1924 &      nbout,tbf_2)
1925      IF (l_dbg) THEN
1926        WRITE(*,*) &
1927 &       "histwrite: 3.4a nbout : ",nbin,nbout,W_F(idf)%W_V(iv)%sopp(io)
1928      ENDIF
1929!-
1930      nbin = nbout
1931      nbout = W_F(idf)%W_V(iv)%datasz_max
1932      CALL mathop(W_F(idf)%W_V(iv)%sopp(io+1),nbin,tbf_2, &
1933 &      missing_val,nbindex,nindex,W_F(idf)%W_V(iv)%scal(io+1), &
1934 &      nbout,tbf_1)
1935      IF (l_dbg) THEN
1936        WRITE(*,*) &
1937 &     "histwrite: 3.4b nbout : ",nbin,nbout,W_F(idf)%W_V(iv)%sopp(io+1)
1938      ENDIF
1939    ENDDO
1940!-
1941!   3.5 Zoom into the data
1942!-
1943    IF (l_dbg) THEN
1944      WRITE(*,*) &
1945 &     "histwrite: 3.5 size(tbf_1) : ",SIZE(tbf_1)
1946      WRITE(*,*) &
1947 &     "histwrite: 3.5 slab in X :", &
1948 &     W_F(idf)%W_V(iv)%zorig(1),W_F(idf)%W_V(iv)%zsize(1)
1949      WRITE(*,*) &
1950 &     "histwrite: 3.5 slab in Y :", &
1951 &     W_F(idf)%W_V(iv)%zorig(2),W_F(idf)%W_V(iv)%zsize(2)
1952      WRITE(*,*) &
1953 &     "histwrite: 3.5 slab in Z :", &
1954 &     W_F(idf)%W_V(iv)%zorig(3),W_F(idf)%W_V(iv)%zsize(3)
1955      WRITE(*,*) &
1956 &     "histwrite: 3.5 slab of input:", &
1957 &     W_F(idf)%W_V(iv)%scsize(1), &
1958 &     W_F(idf)%W_V(iv)%scsize(2), &
1959 &     W_F(idf)%W_V(iv)%scsize(3)
1960    ENDIF
1961!---
1962!-- We have to consider blocks of contiguous data
1963!---
1964    nx=MAX(W_F(idf)%W_V(iv)%zsize(1),1)
1965    ny=MAX(W_F(idf)%W_V(iv)%zsize(2),1)
1966    nz=MAX(W_F(idf)%W_V(iv)%zsize(3),1)
1967    IF     (     (W_F(idf)%W_V(iv)%zorig(1) == 1) &
1968 &          .AND.(   W_F(idf)%W_V(iv)%zsize(1) &
1969 &                == W_F(idf)%W_V(iv)%scsize(1)) &
1970 &          .AND.(W_F(idf)%W_V(iv)%zorig(2) == 1) &
1971 &          .AND.(   W_F(idf)%W_V(iv)%zsize(2) &
1972 &                == W_F(idf)%W_V(iv)%scsize(2))) THEN
1973      kt = (W_F(idf)%W_V(iv)%zorig(3)-1)*nx*ny
1974      tbf_2(1:nx*ny*nz) = tbf_1(kt+1:kt+nx*ny*nz)
1975    ELSEIF (     (W_F(idf)%W_V(iv)%zorig(1) == 1) &
1976 &          .AND.(   W_F(idf)%W_V(iv)%zsize(1) &
1977 &                == W_F(idf)%W_V(iv)%scsize(1))) THEN
1978      kc = -nx*ny
1979      DO kz=W_F(idf)%W_V(iv)%zorig(3),W_F(idf)%W_V(iv)%zorig(3)+nz-1
1980        kc = kc+nx*ny
1981        kt = ( (kz-1)*W_F(idf)%W_V(iv)%scsize(2) &
1982 &            +W_F(idf)%W_V(iv)%zorig(2)-1)*nx
1983        tbf_2(kc+1:kc+nx*ny) = tbf_1(kt+1:kt+nx*ny)
1984      ENDDO
1985    ELSE
1986      kc = -nx
1987      DO kz=W_F(idf)%W_V(iv)%zorig(3),W_F(idf)%W_V(iv)%zorig(3)+nz-1
1988        DO ky=W_F(idf)%W_V(iv)%zorig(2),W_F(idf)%W_V(iv)%zorig(2)+ny-1
1989          kc = kc+nx
1990          kt = ((kz-1)*W_F(idf)%W_V(iv)%scsize(2)+ky-1) &
1991 &            *W_F(idf)%W_V(iv)%scsize(1) &
1992 &            +W_F(idf)%W_V(iv)%zorig(1)-1
1993          tbf_2(kc+1:kc+nx) = tbf_1(kt+1:kt+nx)
1994        ENDDO
1995      ENDDO
1996    ENDIF
1997!-
1998!-- 4.0 Get the min and max of the field
1999!-
2000    IF (l_dbg) THEN
2001      WRITE(*,*) "histwrite: 4.0 tbf_1",idf,iv, &
2002 &      TRIM(tmp_opp),' ---- ',LEN_TRIM(tmp_opp),nbindex
2003    ENDIF
2004!-
2005    IF (W_F(idf)%W_V(iv)%hist_calc_rng) THEN
2006      W_F(idf)%W_V(iv)%hist_minmax(1) = &
2007 &      MIN(W_F(idf)%W_V(iv)%hist_minmax(1), &
2008 &      MINVAL(tbf_2(1:tsz),MASK=tbf_2(1:tsz) /= missing_val))
2009      W_F(idf)%W_V(iv)%hist_minmax(2) = &
2010 &      MAX(W_F(idf)%W_V(iv)%hist_minmax(2), &
2011 &      MAXVAL(tbf_2(1:tsz),MASK=tbf_2(1:tsz) /= missing_val))
2012    ENDIF
2013!-
2014!-- 5.0 Do the operations if needed. In the case of instantaneous
2015!--     output we do not transfer to the time_buffer.
2016!-
2017    IF (l_dbg) THEN
2018      WRITE(*,*) "histwrite: 5.0 idf : ",idf," iv : ",iv," tsz : ",tsz
2019    ENDIF
2020!-
2021    IF (     (TRIM(tmp_opp) /= "inst") &
2022 &      .AND.(TRIM(tmp_opp) /= "once") ) THEN
2023      CALL moycum(tmp_opp,tsz,W_F(idf)%W_V(iv)%t_bf, &
2024 &           tbf_2,W_F(idf)%W_V(iv)%nb_opp)
2025    ENDIF
2026!-
2027    W_F(idf)%W_V(iv)%last_opp = pitau
2028    W_F(idf)%W_V(iv)%nb_opp = W_F(idf)%W_V(iv)%nb_opp+1
2029!-
2030  ENDIF
2031!-
2032! 6.0 Write to file if needed
2033!-
2034  IF (l_dbg) WRITE(*,*) "histwrite: 6.0",idf
2035!-
2036  IF (do_write) THEN
2037!-
2038    nfid = W_F(idf)%ncfid
2039    nvid = W_F(idf)%W_V(iv)%ncvid
2040!-
2041!-- 6.1 Do the operations that are needed before writting
2042!-
2043    IF (l_dbg) WRITE(*,*) "histwrite: 6.1",idf
2044!-
2045    IF (     (TRIM(tmp_opp) /= "inst") &
2046 &      .AND.(TRIM(tmp_opp) /= "once") ) THEN
2047      t_bnd(1:2) = (/ W_F(idf)%W_V(iv)%last_wrt*W_F(idf)%deltat,rtime /)
2048      rtime = (t_bnd(1)+t_bnd(2))/2.0
2049    ENDIF
2050!-
2051!-- 6.2 Add a value to the time axis of this variable if needed
2052!-
2053    IF (     (TRIM(tmp_opp) /= "l_max") &
2054 &      .AND.(TRIM(tmp_opp) /= "l_min") &
2055 &      .AND.(TRIM(tmp_opp) /= "once") ) THEN
2056!-
2057      IF (l_dbg) WRITE(*,*) "histwrite: 6.2",idf
2058!-
2059      itax  = W_F(idf)%W_V(iv)%t_axid
2060      itime = W_F(idf)%W_V(iv)%nb_wrt+1
2061!-
2062      IF (W_F(idf)%W_V(itax)%tax_last < itime) THEN
2063        iret = NF90_PUT_VAR (nfid,W_F(idf)%W_V(itax)%tdimid, &
2064 &               (/ rtime /),start=(/ itime /),count=(/ 1 /))
2065        IF (W_F(idf)%W_V(itax)%tbndid > 0) THEN
2066          iret = NF90_PUT_VAR (nfid,W_F(idf)%W_V(itax)%tbndid, &
2067 &                 t_bnd,start=(/ 1,itime /),count=(/ 2,1 /))
2068        ENDIF
2069        W_F(idf)%W_V(itax)%tax_last = itime
2070      ENDIF
2071    ELSE
2072      itime=1
2073    ENDIF
2074!-
2075!-- 6.3 Write the data. Only in the case of instantaneous output
2076!       we do not write the buffer.
2077!-
2078    IF (l_dbg) THEN
2079      WRITE(*,*) "histwrite: 6.3",idf,nfid,nvid,iv,itime
2080    ENDIF
2081!-
2082    IF (W_F(idf)%W_V(iv)%scsize(3) == 1) THEN
2083      IF (W_F(idf)%regular) THEN
2084        corner(1:4) = (/ 1,1,itime,0 /)
2085        edges(1:4) = (/ W_F(idf)%W_V(iv)%zsize(1), &
2086 &                      W_F(idf)%W_V(iv)%zsize(2),1,0 /)
2087      ELSE
2088        corner(1:4) = (/ 1,itime,0,0 /)
2089        edges(1:4) = (/ W_F(idf)%W_V(iv)%zsize(1),1,0,0 /)
2090      ENDIF
2091    ELSE
2092      IF (W_F(idf)%regular) THEN
2093        corner(1:4) = (/ 1,1,1,itime /)
2094        edges(1:4) = (/ W_F(idf)%W_V(iv)%zsize(1), &
2095 &                      W_F(idf)%W_V(iv)%zsize(2), &
2096 &                      W_F(idf)%W_V(iv)%zsize(3),1 /)
2097      ELSE
2098        corner(1:4) = (/ 1,1,itime,0 /)
2099        edges(1:4) = (/ W_F(idf)%W_V(iv)%zsize(1), &
2100 &                      W_F(idf)%W_V(iv)%zsize(3),1,0 /)
2101      ENDIF
2102    ENDIF
2103!-
2104    IF (     (TRIM(tmp_opp) /= "inst") &
2105 &      .AND.(TRIM(tmp_opp) /= "once") ) THEN
2106      iret = NF90_PUT_VAR (nfid,nvid,W_F(idf)%W_V(iv)%t_bf, &
2107 &                         start=corner(1:4),count=edges(1:4))
2108    ELSE
2109      iret = NF90_PUT_VAR (nfid,nvid,tbf_2, &
2110 &                         start=corner(1:4),count=edges(1:4))
2111    ENDIF
2112!-
2113    W_F(idf)%W_V(iv)%last_wrt = pitau
2114    W_F(idf)%W_V(iv)%nb_wrt = W_F(idf)%W_V(iv)%nb_wrt+1
2115    W_F(idf)%W_V(iv)%nb_opp = 0
2116!---
2117!   After the write the file can be synchronized so that no data is
2118!   lost in case of a crash. This feature gives up on the benefits of
2119!   buffering and should only be used in debuging mode. A flag is
2120!   needed here to switch to this mode.
2121!---
2122!   iret = NF90_SYNC (nfid)
2123!-
2124  ENDIF
2125!----------------------------
2126END SUBROUTINE histwrite_real
2127!===
2128SUBROUTINE histvar_seq (idf,pvarname,idv)
2129!---------------------------------------------------------------------
2130!- This subroutine optimize the search for the variable in the table.
2131!- In a first phase it will learn the succession of the variables
2132!- called and then it will use the table to guess what comes next.
2133!- It is the best solution to avoid lengthy searches through array
2134!- vectors.
2135!-
2136!- ARGUMENTS :
2137!-
2138!- idf      : id of the file on which we work
2139!- pvarname : The name of the variable we are looking for
2140!- idv      : The var id we found
2141!---------------------------------------------------------------------
2142  IMPLICIT NONE
2143!-
2144  INTEGER,INTENT(in)  :: idf
2145  CHARACTER(LEN=*),INTENT(IN) :: pvarname
2146  INTEGER,INTENT(out) :: idv
2147!-
2148  LOGICAL,SAVE :: learning(nb_files_max)=.TRUE.
2149  INTEGER,SAVE :: overlap(nb_files_max) = -1
2150  INTEGER,SAVE :: varseq(nb_files_max,nb_var_max*3)
2151  INTEGER,SAVE :: varseq_len(nb_files_max) = 0
2152  INTEGER,SAVE :: varseq_pos(nb_files_max)
2153  INTEGER,SAVE :: varseq_err(nb_files_max) = 0
2154  INTEGER      :: ib,sp,nn,pos
2155  CHARACTER(LEN=70) :: str70
2156  LOGICAL :: l_dbg
2157!---------------------------------------------------------------------
2158  CALL ipsldbg (old_status=l_dbg)
2159!-
2160  IF (l_dbg) THEN
2161    WRITE(*,*) 'histvar_seq, start of the subroutine :',learning(idf)
2162  ENDIF
2163!-
2164  IF (learning(idf)) THEN
2165!-
2166!-- 1.0 We compute the length over which we are going
2167!--     to check the overlap
2168!-
2169    IF (overlap(idf) <= 0) THEN
2170      IF (W_F(idf)%n_var > 6) THEN
2171        overlap(idf) = W_F(idf)%n_var/3*2
2172      ELSE
2173        overlap(idf) = W_F(idf)%n_var
2174      ENDIF
2175    ENDIF
2176!-
2177!-- 1.1 Find the position of this string
2178!-
2179    CALL find_str (W_F(idf)%W_V(1:W_F(idf)%n_var)%v_name,pvarname,pos)
2180    IF (pos > 0) THEN
2181      idv = pos
2182    ELSE
2183      CALL ipslerr (3,"histvar_seq", &
2184 &      'The name of the variable you gave has not been declared', &
2185 &      'You should use subroutine histdef for declaring variable', &
2186 &      TRIM(pvarname))
2187    ENDIF
2188!-
2189!-- 1.2 If we have not given up we store the position
2190!--     in the sequence of calls
2191!-
2192    IF (varseq_err(idf) >= 0) THEN
2193      sp = varseq_len(idf)+1
2194      IF (sp <= nb_var_max*3) THEN
2195        varseq(idf,sp) = idv
2196        varseq_len(idf) = sp
2197      ELSE
2198        CALL ipslerr (2,"histvar_seq",&
2199 &       'The learning process has failed and we give up. '// &
2200 &       'Either you sequence is',&
2201 &       'too complex or I am too dumb. '// &
2202 &       'This will only affect the efficiency',&
2203 &       'of your code. Thus if you wish to save time'// &
2204 &       ' contact the IOIPSL team. ')
2205        WRITE(*,*) 'The sequence we have found up to now :'
2206        WRITE(*,*) varseq(idf,1:sp-1)
2207        varseq_err(idf) = -1
2208      ENDIF
2209!-
2210!---- 1.3 Check if we have found the right overlap
2211!-
2212      IF (varseq_len(idf) >= overlap(idf)*2) THEN
2213!-
2214!------ We skip a few variables if needed as they could come
2215!------ from the initialisation of the model.
2216!-
2217        DO ib = 0,sp-overlap(idf)*2
2218          IF ( learning(idf) .AND.&
2219            & SUM(ABS(varseq(idf,ib+1:ib+overlap(idf)) -&
2220            & varseq(idf,sp-overlap(idf)+1:sp))) == 0 ) THEN
2221            learning(idf) = .FALSE.
2222            varseq_len(idf) = sp-overlap(idf)-ib
2223            varseq_pos(idf) = overlap(idf)+ib
2224            varseq(idf,1:varseq_len(idf)) = &
2225 &            varseq(idf,ib+1:ib+varseq_len(idf))
2226          ENDIF
2227        ENDDO
2228      ENDIF
2229    ENDIF
2230  ELSE
2231!-
2232!-- 2.0 Now we know how the calls to histwrite are sequenced
2233!--     and we can get a guess at the var ID
2234!-
2235    nn = varseq_pos(idf)+1
2236    IF (nn > varseq_len(idf)) nn = 1
2237!-
2238    idv = varseq(idf,nn)
2239!-
2240    IF (TRIM(W_F(idf)%W_V(idv)%v_name) /= TRIM(pvarname)) THEN
2241      CALL find_str (W_F(idf)%W_V(1:W_F(idf)%n_var)%v_name,pvarname,pos)
2242      IF (pos > 0) THEN
2243        idv = pos
2244      ELSE
2245        CALL ipslerr (3,"histvar_seq", &
2246 &  'The name of the variable you gave has not been declared',&
2247 &  'You should use subroutine histdef for declaring variable', &
2248 &  TRIM(pvarname))
2249      ENDIF
2250      varseq_err(idf) = varseq_err(idf)+1
2251    ELSE
2252!-
2253!---- We only keep the new position if we have found the variable
2254!---- this way. This way an out of sequence call to histwrite does
2255!---- not defeat the process.
2256!-
2257      varseq_pos(idf) = nn
2258    ENDIF
2259!-
2260    IF (varseq_err(idf) >= 10) THEN
2261      WRITE(str70,'("for file ",I3)') idf
2262      CALL ipslerr (2,"histvar_seq", &
2263 &  'There were 10 errors in the learned sequence of variables',&
2264 &  str70,'This looks like a bug, please report it.')
2265         varseq_err(idf) = 0
2266    ENDIF
2267  ENDIF
2268!-
2269  IF (l_dbg) THEN
2270    WRITE(*,*) &
2271 &   'histvar_seq, end of the subroutine :',TRIM(pvarname),idv
2272  ENDIF
2273!-------------------------
2274END SUBROUTINE histvar_seq
2275!===
2276SUBROUTINE histsync (idf)
2277!---------------------------------------------------------------------
2278!- This subroutine will synchronise all
2279!- (or one if defined) opened files.
2280!-
2281!- VERSION
2282!-
2283!---------------------------------------------------------------------
2284  IMPLICIT NONE
2285!-
2286! idf  : optional argument for fileid
2287  INTEGER,INTENT(in),OPTIONAL :: idf
2288!-
2289  INTEGER :: ifile,iret,i_s,i_e
2290!-
2291  LOGICAL :: l_dbg
2292!---------------------------------------------------------------------
2293  CALL ipsldbg (old_status=l_dbg)
2294!-
2295  IF (l_dbg) THEN
2296    WRITE(*,*) "->histsync"
2297  ENDIF
2298!-
2299  IF (PRESENT(idf)) THEN
2300    IF ( (idf >= 1).AND.(idf <= nb_files_max) ) THEN
2301      IF (W_F(idf)%ncfid > 0) THEN
2302        i_s = idf
2303        i_e = idf
2304      ELSE
2305        i_s = 1
2306        i_e = 0
2307        CALL ipslerr (2,'histsync', &
2308 &       'Unable to synchronise the file :','probably','not opened')
2309      ENDIF
2310    ELSE
2311      CALL ipslerr (3,'histsync','Invalid file identifier',' ',' ')
2312    ENDIF
2313  ELSE
2314    i_s = 1
2315    i_e = nb_files_max
2316  ENDIF
2317!-
2318  DO ifile=i_s,i_e
2319    IF (W_F(ifile)%ncfid > 0) THEN
2320      IF (l_dbg) THEN
2321        WRITE(*,*) '  histsync - synchronising file number ',ifile
2322      ENDIF
2323      iret = NF90_SYNC(W_F(ifile)%ncfid)
2324    ENDIF
2325  ENDDO
2326!-
2327  IF (l_dbg) THEN
2328    WRITE(*,*) "<-histsync"
2329  ENDIF
2330!----------------------
2331END SUBROUTINE histsync
2332!===
2333SUBROUTINE histclo (idf)
2334!---------------------------------------------------------------------
2335!- This subroutine will close all (or one if defined) opened files
2336!-
2337!- VERSION
2338!-
2339!---------------------------------------------------------------------
2340  IMPLICIT NONE
2341!-
2342! idf  : optional argument for fileid
2343  INTEGER,INTENT(in),OPTIONAL :: idf
2344!-
2345  INTEGER :: ifile,nfid,nvid,iret,iv,i_s,i_e
2346  LOGICAL :: l_dbg
2347!---------------------------------------------------------------------
2348  CALL ipsldbg (old_status=l_dbg)
2349!-
2350  IF (l_dbg) THEN
2351    WRITE(*,*) "->histclo"
2352  ENDIF
2353!-
2354  IF (PRESENT(idf)) THEN
2355    IF ( (idf >= 1).AND.(idf <= nb_files_max) ) THEN
2356      IF (W_F(idf)%ncfid > 0) THEN
2357        i_s = idf
2358        i_e = idf
2359      ELSE
2360        i_s = 1
2361        i_e = 0
2362        CALL ipslerr (2,'histclo', &
2363 &       'Unable to close the file :','probably','not opened')
2364      ENDIF
2365    ELSE
2366      CALL ipslerr (3,'histclo','Invalid file identifier',' ',' ')
2367    ENDIF
2368  ELSE
2369    i_s = 1
2370    i_e = nb_files_max
2371  ENDIF
2372!-
2373  DO ifile=i_s,i_e
2374    IF (W_F(ifile)%ncfid > 0) THEN
2375      IF (l_dbg) THEN
2376        WRITE(*,*) '  histclo - closing specified file number :',ifile
2377      ENDIF
2378      nfid = W_F(ifile)%ncfid
2379      iret = NF90_REDEF(nfid)
2380!-----
2381!---- 1. Loop on the number of variables to add some final information
2382!-----
2383      IF (l_dbg) THEN
2384        WRITE(*,*) '  Entering loop on vars : ',W_F(ifile)%n_var
2385      ENDIF
2386      DO iv=1,W_F(ifile)%n_var
2387!------ Extrema
2388        IF (W_F(ifile)%W_V(iv)%hist_wrt_rng) THEN
2389          IF (l_dbg) THEN
2390            WRITE(*,*) 'min value for file :',ifile,' var n. :',iv, &
2391 &                     ' is : ',W_F(ifile)%W_V(iv)%hist_minmax(1)
2392            WRITE(*,*) 'max value for file :',ifile,' var n. :',iv, &
2393 &                     ' is : ',W_F(ifile)%W_V(iv)%hist_minmax(2)
2394          ENDIF
2395          IF (W_F(ifile)%W_V(iv)%hist_calc_rng) THEN
2396!---------- Put the min and max values on the file
2397            nvid = W_F(ifile)%W_V(iv)%ncvid
2398            IF (W_F(ifile)%W_V(iv)%v_typ == hist_r8) THEN
2399              iret = NF90_PUT_ATT(nfid,nvid,'valid_min', &
2400 &                     REAL(W_F(ifile)%W_V(iv)%hist_minmax(1),KIND=8))
2401              iret = NF90_PUT_ATT(nfid,nvid,'valid_max', &
2402 &                     REAL(W_F(ifile)%W_V(iv)%hist_minmax(2),KIND=8))
2403            ELSE
2404              iret = NF90_PUT_ATT(nfid,nvid,'valid_min', &
2405 &                     REAL(W_F(ifile)%W_V(iv)%hist_minmax(1),KIND=4))
2406              iret = NF90_PUT_ATT(nfid,nvid,'valid_max', &
2407 &                     REAL(W_F(ifile)%W_V(iv)%hist_minmax(2),KIND=4))
2408            ENDIF
2409          ENDIF
2410        ENDIF
2411!------ Time-Buffers
2412        IF (ASSOCIATED(W_F(ifile)%W_V(iv)%t_bf)) THEN
2413          DEALLOCATE(W_F(ifile)%W_V(iv)%t_bf)
2414        ENDIF
2415!------ Reinitialize the sizes
2416        W_F(ifile)%W_V(iv)%datasz_in(:) = -1
2417        W_F(ifile)%W_V(iv)%datasz_max = -1
2418      ENDDO
2419!-----
2420!---- 2. Close the file
2421!-----
2422      IF (l_dbg) WRITE(*,*) '  close file :',nfid
2423      iret = NF90_CLOSE(nfid)
2424      W_F(ifile)%ncfid = -1
2425      W_F(ifile)%dom_id_svg = -1
2426    ENDIF
2427  ENDDO
2428!-
2429  IF (l_dbg) THEN
2430    WRITE(*,*) "<-histclo"
2431  ENDIF
2432!---------------------
2433END SUBROUTINE histclo
2434!===
2435SUBROUTINE ioconf_modname (str)
2436!---------------------------------------------------------------------
2437!- This subroutine allows to configure the name
2438!- of the model written into the file
2439!---------------------------------------------------------------------
2440  IMPLICIT NONE
2441!-
2442  CHARACTER(LEN=*),INTENT(IN) :: str
2443!---------------------------------------------------------------------
2444  IF (.NOT.lock_modname) THEN
2445    model_name = str(1:MIN(LEN_TRIM(str),80))
2446    lock_modname = .TRUE.
2447  ELSE
2448    CALL ipslerr (2,"ioconf_modname", &
2449   &  'The model name can only be changed once and only', &
2450   &  'before it is used. It is now set to :',model_name)
2451  ENDIF
2452!----------------------------
2453END SUBROUTINE ioconf_modname
2454!-
2455!===
2456!-
2457!-----------------
2458END MODULE histcom
Note: See TracBrowser for help on using the repository browser.