Context Navigation

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

Last change on this file since 1519 was 1517, checked in by mmaipsl, 13 years ago
Add new function histglobal_attr to give GLOBAL ATTRIBUTES in history files. Use this function before histend call.
Property svn:keywords set to `Id`
File size: 80.0 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: