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

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

JB: Fix a bug introduced in 2003 (!)

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