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

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

New handling for time axes

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